At present, most oilfields in eastern China have entered the high water cut stage, traditional sedimentary microfacies characterization cannot meet the production requirements, and fine characterization of single sand body and its internal architecture has become the hot spot in reservoir research^[1,2]. Reserves in fluvial facies reservoirs account for more than 40% of the total reserves of developed oilfields in eastern China^[3]. As an important type of fluvial reservoir, the characterization of meandering river single sand body and its internal architecture has important guiding meaning for future study of remaining oil at the late period of high water cut stage, So the predecessors have done a lot of research and obtained relatively perfect single sand body and architecture analysis methods of meandering river point bar^[4,5,6]. However, due to the complexity of sand cutting and overlapping, the analysis results of meandering point bars are often multi-solution. There are still some uncertainties about how to identify the direction of abandoned channels and predict the internal architectures of point bars accurately^[7]. In this paper, the single sand body and its internal architecture of meandering river point bar are quantitatively characterized by synthesizing the techniques of abandoned channel trend judging, sand body plane combination and empirical formula analysis with horizontal and directional well data. On this basis, a three-dimensional prototype model of three-level architecture interface is established, and the spatial distributions of architecture and remaining oil in meandering river are also studied.

The study area is located in the western part of Shijiutuo bulge of Chengning uplift in Bohai Bay Basin, and is a large- scale drape structure developed in the background of pre-Tertiary paleo-uplift complicated by faults. The oil-bearing strata are mainly in the lower part of Minghuazhen Formation (Nm) and the upper part of Guantao Formation (Ng) of Neogene. The target zone is located in the second sublayer of the third layers group in the lower part of Minghuazhen Formation (NmⅢ₂), which was deposited in meandering river sedimentary environment, and has strip-like channels and sand layers of 6-12 m thick each (Fig. 1). Because of the shallow burial depth, the reservoir is under compacted and weak in diagenesis, and thus having good physical properties, with an average porosity of 35% and an average permeability of 1750×10^-3 μm². Currently, the study area has entered the stage of high water cut development, with local water cut of over 90%, but the recovery rate is only 21%, with a lot of movable oil remaining in the reservoir still. Previous studies have shown that in the development of high water cut stage, the interlayer in reservoir has become an important factor affecting fluid movement and remaining oil distribution^[5,6]. Therefore, the fine anatomy of structure of the sand body and the establishment of related point bar architecture model are helpful to better understand the reservoir heterogeneity and the distribution of remaining oil underground.

The internal architecture anatomy of meandering river reservoirs can be divided into two levels: (1) identifying the sand bodies of each origin inside the meandering river; and (2) studying the occurrence and scale of three-levels architecture elements in point bar^[5,8]. The former is the basis of architectural anatomy, including the identification of point bars and the characterization of abandoned channels. Especially, the description of abandoned channels is directly related to the judgement of the tendency of lateral accretion beddings, which needs to be focused on.

Because abandoned channel is the boundary of the last lateral migration of point bar, it forms an "encirclement" shape with the point bar. Therefore, the method of “defining boundary by abandoned channel and point bar by thickness” is used to describe the plane distribution of the single sand body inside the meandering river^[9].

Point bar is the basic unit and the main body of the "dual structure" of the meandering river, and the vertical positive rhythm of single well is its main feature. Spontaneous potential (SP) curves are bell-shaped or box-shaped, and the amplitudes of resistivity curves vary greatly.

The abandoned channel usually develops around point bars, and thus is an important identification mark of the point bar boundary. According to the different ways of abandonment, abandoned channels can be divided into gradual abandonment type and sudden abandonment type. The gradual abandoned type reflects that the hydrodynamic power gradually weakened. On log curve, they gradually decrease from the high amplitude at the bottom to the small dentate shape. The sudden abandonment type is formed in the sedimentary environment where hydrodynamic force suddenly weakens, and on log curve, it changes rapidly from the high amplitude at the bottom to the small teeth shape near the mudstone baseline^[10,11,12]. Abandoned channels in the study area are mainly sudden abandonment type, reflecting the abrupt change of hydrodynamic sedimentary environment.

The formation process of point bar is an obvious process of "concave erosion and convex accumulation". The sand body of point bar is lenticular in the composite channel and separated from the abandoned channel around it. Therefore, the distribution of point bars and abandoned channels can be effectively identified by sandstone thickness map. Five point bars can be identified in the study area, which are arranged from east to west. The abandoned channel occurring around the point bar is an important identification mark for the spacing of different point bars (Fig. 2).

Plane combination of point bar and abandoned channel is the key to analyze the distribution of single sand body in meandering river. At present, the mainstream method is to combine the distribution pattern of meandering river sand body with dense well pattern. However, this method requires a high density of wells. Because the width of abandoned channel in traditional well pattern is often close to or less than well spacing, it is difficult to tell the direction and bending amplitude of abandoned channel by traditional well pattern^[9]. Although some researchers put forward "envelope of early-final flow lines" method to assist judgment, the final planar prediction results still have multiple solutions because of the complex evolution of meandering river. For example, there are two different explanatory schemes for the trend of abandoned channels in the study area (Fig. 2). It is not clear which one represents the real flow direction, this problem has restricted the further development of meandering rivers sandbodies anatomy for a long time^[7].

Compared with directional wells, horizontal wells can reveal the lateral changes of reservoirs more accurately. With the development of geosteering technology, lithologic boundaries within several meters above or below the horizontal well trajectory can be effectively identified, which provides a new judgment basis for the combination relationship between sand bodies and abandoned channels^{[13,14,15,16,17]}. For example, the geosteering boundary interpretation result of horizontal well A1 (Fig. 3a) shows that there is a mudstone between two sand bodies. By comparing the pinchout shape on both sides of the sand and the half-moon shape of mudstone river channel, combined with the sedimentary characteristics of the study area, this mudstone is judged to be an abandoned channel held between two point bars. Because the left boundary of the abandoned channel has a bigger inclination angle than the right one, according to the geological form of abandoned channel, the left boundary of the abandoned channel was identified as its concave bank, and the right bank was convex bank.

The geosteering interpretation result of horizontal well, intuitive and reliable, is the direct evidence for judging the trend of abandoned channels. But this method is expensive and has little data available. Therefore, a method calculating the dip angle of sand-mudstone interface by using horizontal well resistivity curve has been proposed, to help determine the trend of abandoned channel, which is a supplement and improvement to the existing sand body-architecture anatomy technology of meandering river.

In the resistivity logging process of horizontal wells, with the movement of the instrument, the logging electrode sonde measures the electrical properties of the reservoir. In homogeneous sandstone formation, the resistivity curve appears as a horizontal line (Fig. 4). When the upper edge of the detection radius of the electrode sonde enters the mudstone (point A), the lithology in the detection radius begins to change, and the resistivity begins to drop. When the lower edge of the detection radius leave the sandstone (point C), the lithology of the formation changes from sand-mudstone to mudstone, and the resistivity stop falling and reaches the lowest value. Therefore, the dip angle of the sand-mudstone interface can be quantitatively characterized by the triangular function of the detection radius, the dip angle and the length of the resistivity falling section. Considering that resistivity electrode sonde will decrease in detection radius with the falling of medium-electric resistivity, so it has different detection radius in the sandstone and mudstone. During resistivity logging, the detection radius in sandstone at point A is expressed as R, that at point C mudstone is r, and the detection radius between point A and point C is a quasi-elliptic structure influenced by the mixing of sandstone and mudstone.

In conclusion, when the well trajectory is horizontal, the dip angle of sand-mudstone interface can be calculated by formula (1). When the well trajectory is not horizontal, formula (1) needs to be corrected according to the combination relationship between the horizontal well trajectory and the sand-mudstone interface (Fig. 5). When the sand-mudstone interface is consistent in tendency with the well trajectory (Fig. 5a-5d), formula (1) should be corrected by formula (2); when the sand-mudstone interface is opposite in tendency with the well trajectory (Fig. 5e-5h), formula (1) need be corrected by formula (3).

The applicable conditions of above method are: (1) The interface of sand-mudstone is abrupt or rapid gradual change in lithology. (2) The strata on both sides of the interface are homogenous. (3) The sand body of the stratum should have a certain thickness, and shouldn’t suddenly pinch out or superimpose in complex manner.

In the study area, there are P16H, P40H and A40H and other types of log resistivity curves, corresponding to the detection radius of 0.3, 0.9, 1.8 m in sandstone (25 Ω•m). Among them, the P16H electrode sonde is too small in detection radius, and greatly affected by the flushing zone. Considering the horizontal wells in the study area are largely drilled along the upper part of the sand body, the detection radius of the A40H is too large for the thickness of the reservoir; Therefore, P40H curves with a moderate detection distance were used primarily to calculate in this study to reduce the possibility that the detection range exceeds the sand body.

3.2.1. Rapid qualitative judgement of concave bank and convex bank of abandoned channel in Well A1

Formula (1) shows that in homogeneous sandstone or mudstone formation, under the limited variation of R and r , the value of θ is negatively correlated with the value of L; that is, the larger the value of L, the smaller the value of θ. Therefore, the concave bank and convex bank of abandoned channel can be identified qualitatively and rapidly by L value.

Taking horizontal well A1 as an example, the resistivity begin to decreases at point A (encountering mudstone), reaches the lowest at point A', keeps straight at point A' to C', then goes up at point C’ (encountering sandstone), and eventually restores straight at point C (Fig. 3b). The research shows that Well A1 encountered abandoned channel in A-C section.

By comparison, the resistivity curve on the left side (A-A') has larger slope than that on the right side (C'-C), and larger L value. It can be qualitatively confirmed that the left side of abandoned channel is larger in occurrence than the right side, and then it can be inferred that the left side is the concave bank and the right side is convex bank. The above result has been verified by the geosteering interpretation result of Well A1 (Fig. 3).

3.2.2. Quantitative determination of concave bank and convex bank of abandoned channel in Well H1

The main body of horizontal well H1 is located in point bar 2, and its tail drilled into the mudstone between point bar 1 and point bar 2 (Fig. 2), which is the favorable development position of abandoned channel. After drilling 180 m sandstone, the horizontal well enters the mudstone. The rapid decreased characteristics of logging curves in this section is suspected to encounter the abandoned channel. The resistivity curve (Fig. 6) reaches the highest value of 25 Ω•m at D point, the lowest value of 3 Ω•m at D' point, the L value is 7.8 m, the R value is 0.9 m, the r value is 0.5 m and the β value is 89.97°. Considering that the horizontal well mainly drills along the upper part of sandstone, it is easier to drill up out of the reservoir, so the value of θ was calculated by formula (2) at 10.2°. As Well H1 is inserted into the mudstone at a large dip on the plane, the real dip angle of the sand-mudstone interface must be much larger, not only much larger than the occurrence of the strata in this area, but also larger than the occurrence of the lateral accretion bedding of 5-6°. In addition, the resistivity curve of Well H1 in sandstone section is relatively straight, with very small variation, which excludes the possibility that the well encounters poor sand bodies such as overbank sand through point bar, and shows that the horizontal well drills directly out from point bar to mudstone. The point bar contact the mudstone at a high dip angle, it should be in lateral contact with the abandoned channel on both sides of the point bar, and the direction of the horizontal well is the direction of lateral accretion beddings. In summary, it is finally concluded that the mudstone section at the tail of horizontal well H1 is the abandoned channel, and the well was finally completed in the abandoned channel. Because the occurrence of the left interface is much larger than that of the lateral accretion bedding, according to the geological shape of abandoned channel, it is judged that the left side is the concave bank and the right side is the convex bank of the abandoned channel. Finally, it is determined that "the Possible Trend 2" in Fig. 2 is the real development direction of the abandoned channel in the study area. Due to the prominent position of point bar 4 relative to point bar 3 and point bar 5, the same conclusion can be drawn by using the "early-final stage river flow lines by combining with sand thickness" method to further verify the above result^[7].

Based on the idea of “defining the boundary by the abandoned channel and location of point bar by sand thickness”, the scope of each point bar was delineated according to the single well analysis of sand thickness and curve characteristics. Then, the abandoned channel was combined along the edge of the point bar according to its distribution pattern. With the profile as verification, the identification and dissection of microfacies in meandering river reservoir of the study area were obtained (Figs. 7 and 8). The results show that there are five point bars arranging in bead string in the study area, with distinctive characteristics of abandoned channel; floodplain occurring around as background facies is mainly flooded mudstone, showing stark distinction from the point bars.

As the most important geomorphic unit of meandering river, point bar is formed by lateral accretion and consists of many lateral accretion beddings and lateral accretion bodies^[18,19,20]. Previous studies also have shown that the point bar is mainly composed of the above two architecture elements. Therefore, quantitative characterization of the architecture parameters of the above architecture elements is the key to fine dissection of point bar^[18].

For the measurement method of lateral accretion beddings occurrence, the research has been relatively mature^{[5-9, 21-22]}. Based on reviewing previous research methods, we used the classical empirical formula method to quantitatively calculate the dip angle of the lateral accretion beddings of point bars in the study area.

A lot of research has been done on the internal architecture pattern of point bars, and a series of empirical formulas have been summarized^{[8-10, 22-24]}. Among them, Leeder^[23] collected 107 meandering river examples and established the formula for the dip angle of the lateral accretion bedding in the point bar under high-bending conditions.

As Leeder’s formula has been widely accepted in the industry, this formula was used to calculate the dip angle of the lateral accretion bedding in the study area. According to the statistical results of abandoned channels, the thickness of the full bank ancient channels is mainly 5.0-6.5 m, so the width of the full bank abandoned channels was estimated at 80-120 m, and the dip angle of the lateral accretion bedding was about 5°.

4.2.1. Empirical formula method

Based on the previous research results^{[6-10, 20]}, the Leeder formula was used to calculate the scale of the lateral accretion bodies in the study area too. In accordance with formula (6), the scale of the lateral accretion body is 2/3 of the width of the full bank channel, and according to the above calculation of full bank width of the abandoned channel, the scale of the lateral accretion bodies in the study area is 55-80 m.

4.2.2. Calculation with horizontal well logging data

By identifying the lateral accretion beddings in horizontal well, the sandstone length between two adjacent lateral accretion beddings can be obtained, that is the apparent width of the lateral accretion body. The true width of the lateral accretion body can be calculated by correcting the apparent width with the angle between the horizontal well and the inclination of the lateral accretion beddings. The horizontal well H2 encountered three lateral accretion beddings and two lateral accretion bodies. The two lateral accretion bodies are 73 m and 75 m wide respectively, and the corresponding angle is 20°. From this, the width of the lateral accretion bodies drilled in the horizontal well were calculated at 68.6 m and 70.5 m respectively (Fig. 9).

In summary, based on empirical formula method and supplemented by horizontal well logging data analysis, the calculation results of the two methods are similar, and the analysis results of horizontal well logging data are within the scope of the calculation results of empirical formula method. Therefore, the width of lateral accretion bodies in the study area is confirmed finally between 55 m and 80 m, and largely near 70 m.

On the basis of the above analysis results, fine architecture dissection was carried out for the five identified point bars. This includes analyzing spatial distribution characteristics of lateral accretion beddings mainly based on the well data with the characteristics of lateral accretion beddings verified by horizontal well, and the scale and occurrence of architecture elements, and mapping the top plane of the lateral accretion beddings in the interior of the point bar by referring to the classical architecture pattern of meandering river^[5] (Fig. 8).

Compared with traditional modeling technology, architecture modeling can simulate the third or fourth level, or even smaller-level architecture interfaces, and greatly improve the accuracy of heterogeneity characterization in reservoir. A lot of work has been done to characterize the three-dimensional architecture of point bars in meandering rivers, many three-dimensional models and corresponding modeling methods have been established to characterize the three-level architecture^{[25,26,27,28,29]}. At present, grids are commonly used to characterize lateral accretion beddings, the models usually have a vast number of grids, making the computation workload huge and time-consuming. If the model is upscaled, the architecture of lateral accretion beddings could miss. These reasons restrict the wide application of lateral accretion beddings model during the oilfield development^[30,31].

In this work, a three-dimensional architecture model was built based on the principle of surface constraint (Fig. 10). Firstly, the related bounding surfaces of point-bar architecture elements (lateral accretion beddings and lateral accretion bodies) were generated based on fine dissection results of reservoir architecture. Then, according to the lateral accretion process of point bar, the top and bottom surfaces of different architecture elements in the point bar were separated, and then they were combined and encapsulated according to their spatial superposition order. Finally, the lateral accretion beddings and the lateral accretion bodies were divided into grids separately, because different lateral accretion beddings and lateral accretion bodies are separated by the relevant surfaces, different mesh accuracy can be designed for different architecture elements in the upscaling process, to solve the contradiction between the big mesh number of the model and the accuracy of the lateral accretion beddings (Fig 11).

Using the above method, the lithofacies model and related physical properties model, which can characterize the third- level architecture in the study area, were established by recombining and overlapping the different point bars in the meandering river (Fig. 12). As a constraint condition of property modeling, lithofacies model is mainly established by deterministic modeling and stochastic modeling method. In the interior of abandoned channel and lateral accretion bedding, deterministic modeling method was adopted to define the dominant lithofacies. For the lateral accretion bodies and flood plain, the sequential indicator simulation method which is suitable for small well spacing and can better simulate complex anisotropy was adopted to effectively characterize the distribution of sandstone in point bars. Finally, taking the lithofacies model as the facies control condition, the related three-dimensional physical property model was established by using the method of sequential Gauss simulation. In the modeling process, the lithology and physical properties of lateral accretion beddings and lateral accretion bodies took the known drilling data strictly, and by setting the relevant variation function, the simulation results of each architecture unit of the final model was ensured to conform to the real underground situation.

All structural units in meandering river are well reproduced in the architecture model (Fig. 12). The meandering river course has obvious tortuosity and point bars are in lateral contact with the abandoned channel. Lateral accretion beddings are developed in the interior of the point bars, and the plane arc characteristic and the inclined profile characteristic of the lateral accretion beddings are also well reproduced. Considering the need to limit the total number of grids (usually no more than 10⁶) in the numerical simulation process, non-uniform upscaling method was adopted to upscale the architecture model. That is, under surfaces constraints, different grid densities were applied to different lateral accretion beddings and bodies, to ensure that the morphological integrity of lateral accretion beddings is preserved while upscaling the sand body model (Fig. 11d). The upscaled model could not only accurately characterize the relevant point bars, abandoned channels and lateral accretion beddings, but also greatly reduce the grid number, with a total effective grids of 78×10⁴ only. The model reached the goal of restoring the distribution pattern of underground architecture with a small number of grids, laying foundation for further reservoir numerical simulation and remaining oil prediction.

The purpose of establishing reservoir architecture model is to predict the distribution of remaining oil. In this work, remaining oil distribution in target sand body was predicted by numerical simulation. The simulation results show that with reasonable number and types of grids, this model has a faster arithmetic speed and better convergence, and is suitable for repeated history matching in oilfield industrial production. Fig. 13 shows predicted distribution of remaining oil by June 2017 from numerical simulation. In order to clearly distinguish different types of architecture units, the grids of abandoned channels and lateral accretion beddings have been hidden.

The simulation results show that at the level of compound meandering river, different point bars are mainly sheltered by abandoned channels, so remaining oil exists where well pattern is not complete. For example, due to the lack of water injection wells, the point bar 4 has higher remaining oil saturation than the other four point bars (Fig. 13a). But within a single point bar, the distribution of remaining oil is more affected by the three-level architecture interface (lateral accretion beddings). The bottom of the channel is mainly a high permeability zone, where the injected water rushes rapidly and the remaining oil saturation is relatively low; because of the blocking effect of lateral accretion beddings, the middle and upper part of the channel are the sites remaining oil exists, and the upper part of the reservoir has higher remaining oil saturation than the middle part (Fig. 13b). Therefore, the upper part of the reservoir in lateral accretion beddings farthest away from the water injection well not drilled by the production well has the most remaining oil.

In this paper, the dip angle of sand-mudstone interface on both sides of abandoned river is calculated qualitatively and quantitatively by using resistivity curve of horizontal well, and then the bending mode of abandoned channel is determined. This is a supplement and improvement to the existing anatomical technology of architecture sand body of meandering river.

Horizontal well H1 encountered abandoned channel deposits at the tail of the horizontal section. Because the occurrence of the left interface is much larger than that of the lateral accretion bedding, according to the geological form of abandoned channel, the left side was identified as the concave bank and the right side the convex bank, and then the true trend of the abandoned channel in the study area was determined.

The width of the lateral accretion bodies of the point bar is mainly 55-80 m, and the lateral accretion beddings have largely dip angles of around 5-6 degrees, and tendency pointing to the development side of the abandoned channel. Extending from the upper to the bottom part of the point bar, they accounts for about 2/3 of the sand thickness of the point bar.

Different from traditional grid constraints, surface restriction and non-uniform upscaling method were used to establish the inner architecture model of point-bars, so the fine 3D construction model of multi-point bar compound channel was built with fewer grids, which allows faster numerical simulation of the whole study area, and serves better for optimizing water injection and remaining oil prediction in the oilfield.

Nomenclature

GR—natural gamma, API;

h—full bank depth of channel, m;

H—vertical distance between two lateral accretion beddings, m;

L—distance between point A and point C, m;

r—detection radius of electrode sonde in mudstone, m;

R—detection radius of electrode sonde in sandstone, m;

R_lld—deep lateral resistivity, Ω•m;

R_lls—shallow lateral resistivity, Ω•m;

R_P40H—horizontal well resistivity, Ω•m;

SP—spontaneous potential, mV;

W—full bank width of channel, m;

W_c—horizontal width of lateral accretion body, m;

α—angle between sand-mudstone interface and well trajectory, (°);

β—inclination of well trajectory, (°);

θ—dip angle of sand-mudstone interface, (°);

γ—dip angle of lateral accretion bedding, (°).