Multiscale Modeling of Meandering Fluvial Reservoir Architecture Based on Multiple-Point Geostatistics: A Case Study of the Minghuazhen Formation, Yangerzhuang Oilfield, Bohai Bay Basin, China

Meandering river reservoirs are essential targets for hydrocarbon exploration, although their characterization can be complex due to their multiscale heterogeneity. Multipoint geostatistics (MPS) has advantages in establishing reservoir architectural models. Training image (TI) stationarity is the main factor limiting the uptake of MPS modeling algorithms in subsurface modeling. A modeling workflow was designed to reproduce the distribution of heterogeneities at different scales in the Miocene Minghuazhen Formation of the Yangerzhuang Oilfield in the Bohai Bay Basin. Two TIs are established for different scales of architecture. An initial unconditional model generated with a process-based simulation method is used as the megascale TI. The mesoscale TI of the lateral accretion layers is characterized by an uneven spatial distribution of mudstone in length, thickness, frequency, and spacing. Models of different scales are combined by the probability cube obtained by lateral accretion azimuthal data as an auxiliary variable. Moreover, the permeability function sets are more suitable than the porosity model for collaboratively simulating the permeability model. Model verification suggests this workflow can accurately realize the multiscale stochastic simulation of channels, point bars, and lateral accretion layers of meandering fluvial reservoirs. The produced model conforms geologically realistically and enables the prediction of interwell permeability variation to enhance oil recovery.


Introduction
Meandering river deposits create significant hydrocarbon reservoirs. However, it is challenging to describe the geometric shapes of architectural elements since these reservoirs are characterized by fast planar phase change, complex sand body stacking, and substantial heterogeneity [1,2]. Describing the architectural elements of these reservoirs has been proven to be an effective procedure for calculating proven reserves and enhancing oil recovery [3]. Megascale reservoir architectures include genetic units and their boundaries, such as the point bar and channels. Its heterogeneity can be found in the genetic unit boundaries, sand body stacking pattern, and connectivity [4][5][6][7][8]. Mesoscale reservoir architectures are zonation of permeability within genetic units and fluid flow baffles within genetic units. Its heterogeneity in fluid flow occurs within point bars, and it is related to the development degree of lateral accretion layers and affects the variation of porosity and permeability [9,10]. Previous studies have focused on the sedimentary structure and lithologic heterogeneity in fluvial depositional systems, but little attention has been focused on the heterogeneity at the architectural element scale, especially the heterogeneity variation in model parameters [11][12][13][14]. Mesoscale heterogeneity in architectural elements becomes critical at extending reservoir production periods due to the increased sensitivity of fluid migration to smaller-scale heterogeneity. Thus, multiscale heterogeneity modeling of meandering river reservoirs is urgently needed to explore and develop meandering river reservoirs. Moreover, reservoir heterogeneity is mainly manifested in the permeability model. The trend of permeability variation of reservoir architecture can be predicted by coring data.
The current modeling approaches have limitations regarding the reproducibility of facies and replicating the input data. The traditional modeling method of sequential indicator simulated (SIS) matches point data (wells) based on only an algorithm without regard to the information on sedimentary process [15][16][17][18]. This method commonly produces unrealistic heterogeneous architectural model elements and overestimates the connectivity of these elements [19]. Object-based modeling (OBM) realistically reproduces heterogeneous architectural details in fluvial facies modeling. However, in an actual model, the connectivity of reservoir elements is easily overconnected by unrealistic shapes [20,21]. Therefore, this method is only suitable for discrete attribute models in simple reservoir elements [22]. Another process-based modeling (PBM) is constrained by physically inspired rules and corresponds to processes occurring in the strata. This method does not allow conditioning on dense well data, but the rules and statistics lifted from such PBM results can be used to parameterize a TI.
Multipoint geostatistics (MPS) creates a new approach to integrate geological concepts for reservoir modeling [23][24][25]. MPS is advantageous for obtaining quantitative deposition information by constructing a training image (TI) containing a geological understanding of the study area. Furthermore, modeling based on MPS captures objects by considering more than two points (cells) simultaneously and reproduces complex morphologies, such as those of the deposits of meandering rivers (e.g., point bars and channel fills) and nonlinear spatial correlations [26]. The design of TIs with the right features at the reservoir grid-scale is essential for MPS facies modeling [21]. The sedimentological features of facies (patterns) should not be confined to specific locations within parts of a TI but should instead be distributed homogeneously and repeated reasonably throughout the entire TI [27,28]. On the other hand, it is challenging to reproduce complex geological trends in the resulting reservoir models, although a complex TI incorporates more data events and more accurate conditional probability because it is nonstationary. Moreover, the larger the size and resolution (number of cells) of the TI, the longer it will take to realize. Therefore, the reproduction of nonstationarity with an evident trend (e.g., direction and compression ratio) is attempted by applying auxiliary variables when performing MPS simulations [29,30]. MPS is suitable for architectural modeling of meandering reservoirs with a relatively specific pattern.
In this paper, a multiscale reservoir architectural modeling workflow is studied, including meandering fluvial channels and point bars with lateral accretion elements. First, the meandering fluvial reservoir architecture was characterized in the Minghuazhen Formation of the Yangerzhuang Oilfield. Second, the obtained geometric architectural parameters that describe the fundamental facies features were used to establish two TIs of different scales. Then, architectural modeling was presented for multiscale reservoirs via the MPS method. The desired trends in facies distribution and orientation were achieved by using auxiliary variables. Finally, the variation function of rock permeability was further used to improve the modeling accuracy.

Geological Setting
The Huanghua Depression is a secondary tectonic unit of the Bohai Bay Basin, situated between the Cangxian uplift and the Chengning uplift (Figure 1(a)). The Yangerzhuang Oilfield studied here is an anticlinal structure located in the Zhaojiapu fracture zone in the Huanghua Depression ( Figure 1(b)). The primary sediment source is from the Chengning uplift in the south. It mainly includes the Neogene Guantao Formation and Minghuazhen Formation ( Figure 2). It is a set of typical reservoirs formed by highcurvature meandering fluvial sedimentary systems in the Minghuazhen Formation. Reservoir types mainly include structural reservoir and lithologic-structural reservoir.
The target formation, NmIII, belongs to the lower part of the Neogene Minghuazhen Formation (Figure 2). This unit comprises weakly diagenetic feldspathic fine sandstone shallowly buried at depths of 1340-1900 m [31]. As the significant oil group of the Yangerzhuang Oilfield, it has a porosity of 26%-34%, a permeability of approximately 2.8 μm 2 , and a water-cut rate generally higher than 94.99% (these data are from oilfield exploration and development sources). The drilled wells are primarily located in the central rift zone, with relatively few wells in surrounding areas. Therefore, integrated analysis of well logs and seismic data is necessary for the relatively sparse well network area due to the uneven distribution of these drilled wells. Previous studies of the Yangerzhuang Oilfield have always focused on the reservoir structure and its effects on the remaining oil, but few quantitative studies have investigated the multiscale heterogeneity in the reservoir [32][33][34].

Data and Methods
3.1. Data Acquisition. Over 9.6 km 2 in the study area, 3D seismic reflectivity data were acquired with a bin spacing of 12:5 m × 12:5 m. These data have a zero-phase polarity. The seismic data have a dominant frequency of approximately 40 Hz sampled at 4 ms, ranging from 8 to 60 Hz, which yields a vertical resolution of approximately 12.5 m.
The exploration wells were drilled during the years 1974-2015. The average well spacing is less than 100 m, while the minimum well spacing is 30 m in the dense well network area. Nine wells were cored extensively during target layers of the reservoir in the dense well network, and over 137 m of conventional cores recovered. Wireline logs from 262 wells were integrated with seismic data to establish 2 Geofluids the stratigraphic framework of the area and architectural boundaries at different scales.

Identification of Reservoir Architectural
Parameters. The 3D seismic reflectivity data and well log data were analyzed in light of the parameters and empirical formulas based on modern sedimentation. The complex reservoir architecture of a meandering river was parameterized on different scales (the megascale and the mesoscale) in the dense well network area. The geometric parameters, such as the thickness and dip angle of the lateral accretion elements, were obtained by the paired wells method, the abandoned channel surface method, and the width-to-depth ratio method. A single channel can be identified by analyzing the depth displacement and the characteristics of the well log curves in channel sand bodies due to differences in the provenance supply, sedimentary environment, and deposition rate. Floodplain or overflow bank deposits were often developed at single-channel boundaries with an abrupt variation in the sand body thickness in the fluvial channel. Then, a   3 Geofluids further division scheme of the single point bar was proposed by analyzing the variation in intercalation thickness, the different characteristics of the well log curves, and the relationships of physical properties in the point bar complex. In addition, dynamic data (such as tracer detection) considering the sedimentary pattern can also be treated as a reference to identify a single point bar.

TIs and Auxiliary Variables Prepared for the MPS Model.
TIs can be obtained from satellite images of modern analogs, human-computer interactive drawings, or unconditional PBM simulation outputs as idealized models representing the fundamental architectural features of the depositional system in the study area. To address the problem of multiscale modeling, the TI should include quantitative information on facies architectures (e.g., facies proportions and channel geometry) at different scales. On the other hand, the TI should be simple enough to reduce the amount of computation. Thus, two TIs of different scales corresponding to the architectures at different scales were developed.
The auxiliary variables were gridded properties that described a spatial variation trend for the modeled property, as the missing geological information was not illustrated in the TIs. They also facilitated reproducing nested architectures, in which simulations were undertaken at multiple scales to predict heterogeneities in a reservoir model. Several auxiliary variables can constrain simulations, such as the distance from the channel middle line, the sand body probability map, and the lateral accretion azimuthal probability cube, and allow the resultant models to account for nonstationarity related to variations in the space of sedimentary features.
It was difficult to detail the sedimentary characteristics of the underground reservoirs in the relatively sparse well network area using only basic well logging information. Thus, it was unreliable to use only the sedimentary facies for modeling constraints. Therefore, the permeability distribution function was introduced as an auxiliary condition and was combined with the sedimentary facies and reservoir stratification to achieve high-accuracy prediction.

Workflow of Establishing a Multiscale Reservoir Model.
A multiscale 3D reservoir model of the meandering fluvial system was constructed with a hierarchical workflow ( Figure 3). First, the 3D TI suitable for the megascale was obtained with the unconditional PBM simulation based on an analysis of the geometric parameters in each sedimentary element [35]. Additionally, the mesoscale 3D TI inside the point bar was drawn interactively based on the characteristics of the lateral accretion elements. Second, the facies map was obtained from the seismic interpretation and was consistent with well logs, which could summarize the lateral accretion azimuthal data. Third, the MPS method was applied in a multiscale architectural model of the meandering fluvial system reservoir with constraints of the seismic attribute cube and the lateral accretion azimuthal probability model. Finally, the property model was cooperatively simulated based on the distribution functions regressed from the core physical property data.  4 Geofluids based on grain size, thickness, and sedimentary structure in the meandering fluvial system (Table 1; Figure 4): one conglomerate facies (Cc), four sandstone facies (Sm, Sp, St, and Shv), two fine-grained siltstone facies (Fr and Fh), and one mudstone facies (Mm). The terminologies are referred from Miall [22]. A vertical lithofacies sequence of an ideal meandering river model from bottom to top is channel lag deposits above an erosion surface, massive bedded sandstone, cross-bedded sandstone, horizontally laminated sandstone, ripple laminated siltstone, and massive mudstone ( Figure 5). The meandering channel migrates laterally due to hydrological variation and forms lateral accretion elements and multiple secondary lithologic cycles in modern sedimentary environments. In addition, abandoned channels develop sporadically in the meandering belt. Thus, this process results in heterogeneous multicycle variations in lithology and physical properties [22,41].

Petrophysical Characteristics of Cores.
The experimental data of the wells illustrate that moderate correlation is observed between the vertical sedimentary cycle and the lateral units. Moreover, the experimental data from the coring wells in different meandering belts, point bars, or sedimentary facies were compared. The wells have almost the same spatial distribution regularities in adjacent meandering belts. A moderate positive correlation is seen between the median grain size and the depth, corresponding to a linear increase ( Figure 6(a)). Similarly, some moderate correlation is observed between the shale content and the depth, corresponding to a linear decrease ( Figure 6(b)). Although it has a similar linear declining trend, the shale content in the point bars is significantly lower than that in the channels. A moderate linear correlation is also seen between the logarithmic values of the permeability and the median grain size ( Figure 6(c)). Finally, a fairly high positive correlation is seen between the permeability and depth, corresponding to an exponential increase ( Figure 6(d)). In addition, the permeability and shale content in point bars increase slightly in the direction of fluvial deposition in a given meandering belt.
The least square method was used to fit the regression functions. The residual analysis was performed with a confidence interval of 75% to optimize the regression functions. As a result, the vertical distribution functions of the median grain size (1), shale content of the point bar (2), shale content of the channel (3), and permeability (4) where x is the distance from the top of the layers and R 2 is the goodness of fit.
The observed distribution in the data may reflect little difference between channel sand and point-bar sand bodies in the exponential function of permeability and median grain size (Figure 6(c)). Furthermore, the relationships likely indicate that the vertical exponential function of permeability is mainly controlled by the thickness of the sand bodies in the fluvial channel ( Figure 6(d)). Two regression lines are formed in the dataset, representing the physical property distribution at the top and bottom of the sand bodies in different periods ( Figure 6). This quantification finds application in subsurface workflows, where it can provide some geological constraints for static models requiring architectures to be realistic.
The porosity shows a noticeable linear increase at the top of a point bar but remains almost unchanged from the middle to lower parts (Figure 7(a)). Nevertheless, the correlations between the porosity and the permeability of the sand bodies in the meandering fluvial channels are relatively modest (Figure 7(b)). Relatively modest correlations are also observed between porosity and the logarithmic values of permeability in point-bar sand bodies, ignoring porosity values less than 25% (Figure 7(b)).
Because the coring well is drilled vertically through several lateral units in the point bar, mudstone and sandstone are interbedded in the cores, as in the sample wells in Figure 8. The permeability distribution function of the lateral accretion elements in the point bar can be fitted after ignoring the mudstone value (Figure 9(a)). The petrophysical variation trend is almost identical under the same geological and hydrological characteristics in adjacent lateral units. As a result, the measurement and analysis data can be regarded as data from different positions in the same lateral accretion element ( Figure 8). The permeability increases approximately exponentially from the top (mudstone value) of the point bar to the bottom (sandstone value) (Figure 9(b)). This might be due to the mudstone covering the top of the lateral accretion elements, while the sand body connects the lower part. This is the location where the lateral accretion layers fade away, approximately 1/2 or 1/3 of the upper segment of the point bar, where the permeability increases rapidly (Figure 9(a)).
The observed distribution in the data may partially reflect the fact that the distribution trends of permeability above and below the lateral accretion layer are different. This might be because the original sediments below the lateral accretion layer were eroded to different degrees by later river stages, resulting in the differential preservation of the sediments above and below the lateral accretion layer [42]. The permeability distribution along the lateral accretion elements is according to the exponential function set: where P is the permeability of the lateral accretion elements, x is the distance to the top of the point bar, and C m is the permeability of mudstone, which is generally less than 20 mD. It can be determined according to the actual minimum permeability value in different point bars. L is the 5 Geofluids Channel lag deposits are coarse-grained bedload deposits in a high-energy fluvial flow [22]. clast material present sporadically [22,37] Planar cross-bedded sandstone (Sp) Fine-grained to medium-grained arenite~1 The low-angle cross-bedding is less than 10°, which means moderate-energy flow. The architectural interface has prominent lithologic and electrical characteristics [33,44]. The gamma ray (GR) log curve features clear returns, the microelectrode log curve has low values, and the compensated neutron log curve features significant increases in values ( Figure 10). Reservoir sedimentary characteristics, such as the sand body shape and combination, development scale, and connectivity, are summarized in the study area [1]. The architectural (facies) boundary range can be identified by well data in the dense well network area. However, this is difficult to determine in the relatively sparse well network area and has to be derived with the empirical formula.

Occurs
Point bars and channels are usually formed by the cutting and overlapping actions of fluvial channels. As a point-bar complex, a single channel is formed by meandering fluvial channel migration (Figure 8). The abandoned channel represents the end of a point bar and indicates the diversion of the fluvial channel. Thus, the thickness of a sand body determines the point bar complex location, while the abandoned channel determines its boundary on the plane. On the other hand, in the vertical direction, a compound rhythm is formed by the multiple positive rhythms of an upward-fining cycle caused by the superposition of point bars [1,45].
Predecessors fitted the curvature radius, length of the point bar, and bankfull width of the river with the least square method and obtained the empirical formula. The curvature radius and length limits of the point bar were determined with a 75% confidence level [46,47]. The architectural parameters of the meandering fluvial reservoir are obtained from the calculation shown in Table 2. The thickness of active or abandoned channel fill and barshaped architectural elements-where these elements are fully preserved-acts as a direct indicator of the maximum bankfull depth (h) with compaction correction (Figure 8) [37,47]. The relationship formula W r = 6:8h 1:54 was proposed by Leeder [37] to establish the relationship between bankfull depth (h) and width (W r ). The width of the meander belt (W m ) and the width of the point bar (W d ), as shown in Table 2, are obtained using the empirical formulas below.

Architecture Elements inside the Point Bar.
The argillaceous, calcareous, and physical intercalations are the primary sources of heterogeneity and are treated as the architecture elements developed inside the point bar. Argillaceous intercalations composed of mudstone or silty mudstone are developed in the study area and mainly formed during the intermittent period of flooding [1]. It is determined by well log characteristics that the return of the microelectrode log curve is over 30%, and the amplitude variations of the micronormal log and the microinverse log are almost zero. The lithology of the calcareous intercalations is dominated by calcareous siltstone with a thickness of approximately 15-20 cm. This is shown as an apparent peak of the microelectrode log and the three porosity logs. Furthermore, physical intercalations with low porosity and permeability are rare in the study area due to inadequate compaction.
A connected-well section perpendicular to the river channel is constructed by linking paired wells less than 50 m apart and located in the same point bar in the dense well network area. Therefore, the extended length, dip angle, and horizontal spacing of lateral accretion layers can be quantitatively described. The connected-well section can be employed to calculate the horizontal spacing and the attitude of the lateral accretion layers. The dip angle of the lateral accretion layers can be calculated indirectly with the abandoned surface (the last surface of the lateral accretion sand body before the meandering channel was abandoned) (Figure 8; Table 2).
However, due to the large well spacing and limited data, the bankfull width, bankfull depth, and dip angle of the lat-eral accretion layer in the relatively sparse well network area are calculated using the empirical formula: where θ is the dip angle of the lateral accretion layer (°), h is the bankfull depth of the river channel (m), and W r is the bankfull width of the river channel (m).

Model Grid
Setting. The accuracy of the plane grid is 10 m × 10 m, and the total plane grid is 523 × 339. The vertical grid size was set at 0.2 m to reflect the longitudinal petrophysical variation in the lateral accretion elements. The same grids were applied to all simulations, but different grids for each TI were created and are idealized and scale independent.

Training Images of Different Scales.
The TI is a 3D architectural product that presents reservoir heterogeneity, incorporating the geometries, spatial distribution, and connectivity of geologic bodies [48]. It is treated as a pattern guide to calculate the correlation between point data and helps determine the interwell prediction results of the model.  8 Geofluids systems with high curvature and lateral accretion ( Figure 11) [20]. Parameter sensitivity analysis can help screen suitable TIs with reasonable tests of the location and frequency of sedimentary events following geological rules [35]. The mesoscale TI incorporates only facies representing point-bar, floodplain, and channel-fill deposits, and the active channels are treated as channel-fill deposits (Figure 12(a)). An unconditional PBM simulation generates the mesoscale TI with the use of geometric parameters ( Table 2). The mesoscale TI drawn interactively shows the uneven distribution of the length, thickness, spacing, and continuity of the lateral accretion layers, predicting petrophysical heterogeneity inside the point bar. In addition, it characterizes the connectivity at the bottom of the point bar and the impermeability of the lateral accretion layers (Figure 12(b)).
Therefore, the TIs of this work use only a relative proportion of the lithofacies and the facies relationship to reproduce some high-order properties rather than univariate statistics [1,13].

Auxiliary Variables and Trend
Model. Probability maps have been used to ensure realistic spatial transitions between facies in simulations, particularly the reproduction of geometries and variations in sedimentological features, such as channel-belt and floodplain facies in relatively sparse well networks. Facies analysis of isochronous stratigraphic slices can yield probability maps, and it also helps predict the horizontal distribution of sediments (as shown in the facies maps in Figure 13(b)) [49]. The isochronous seismic events corresponding to the floodplain sediments were selected as relatively reliable isochronal correlation indicators and were linearly interpolated and sampled to extract isochronous stratigraphic slices (Figure 13(a)) [50][51][52]. The amplitude value of sandstone is generally higher than that of mudstone after 90°phase transformation. Therefore, it is simple to interpret high amplitude values as sandstone and low and negative amplitude values as mudstone (Figure 13(a)). However, some amplitude values of sandstone and mudstone overlap.

Geofluids
The lateral accretion azimuthal data control the distribution trend of the lateral accretion layers (Figure 8). Furthermore, since the correlation between porosity and permeability is relatively modest (Figure 7(b)), function (5) is employed to modify the permeability model of the channels. Furthermore, the distribution function set of permeability in the lateral accretion direction constrains the permeability models inside the point bars.

Sedimentary Facies
Modeling. Sedimentary facies are often an essential constraint of multiscale reservoir heterogeneity. Therefore, a hierarchical approach to facies modeling has been taken to employ TIs at differing scales to model the stratal architecture of a particular fluvial meandering system. The expected or desired geological features are forcibly incorporated into the models by using auxiliary variables.
The TI (Figure 12(a)) and probability maps are used to simulate the megascale model with the constraint of sedimentary facies. In contrast, the lateral accretion layers are regarded as facies on the same scale as the point bar when simulating the mesoscale model. Therefore, in this model, the lateral accretion layers in different point bars are com-bined by simply distinguishing the permeable facies and impermeable lateral accretion layers. The lateral accretion azimuthal probability model generated from the lateral accretion azimuthal data (Figure 14(a)), which were normalized to 0 to 360 degrees, is used to constrain the orientation of the lateral accretion layers. According to the abandoned channel and point bar shape, the lateral accretion azimuthal probability model considers individual characteristics of the lateral accretion layers in different point bars. Finally, the mesoscale MPS model is simulated with the TI for lateral accretion layers (Figure 12(b)). Because the distribution position indicates the formation period and shape of the point bar, the abandoned channel becomes a critical connection to match the models at two scales. The lateral accretion layer model at the mesoscale is a stochastic simulation based on the entire megascale model (Figure 14(b)).

Porosity and Permeability
Modeling. Sequential Gaussian simulation (SGS) was used to simulate the distribution of crucial petrophysical parameters such as effective porosity and permeability under the constraint of facies. Each facies corresponds to distinct parameters. For example, the major  10 Geofluids variogram direction of the channel sand bodies is parallel to the fluvial direction, while the minor variogram direction is perpendicular to it. In contrast, the major variogram direction of the point-bar sand bodies is the lateral accretion direction, while the minor variogram direction is perpendicular to it. The details of the variogram adjustment are not described in this paper. These are high-fidelity models based on the available data and eliminate outliers to ensure consistency between petrophysics and electrofacies [53]. The porosity of reservoir sandstone in the target layer is mostly between 19.3% and 37.6%, with an average value of 32%. The porosity distribution shows a linear downward trend at the top of the point bar, while the porosity barely changes from the middle to the lower part. The inflection point of the porosity regression curve is approximately 1/5 down from the top of the sandstone (Figure 7(a)). The porosity range is 20%-38% in the simulated realization, with an average value of 31% (Figure 14(c)), similar to the value range obtained from coring well analysis. Furthermore, the effective porosity of mudstone is much lower than that of sandstone and is less than 5%.
The correlation between porosity and permeability in the target reservoir is relatively modest (Figure 7). The permeability exhibits a log-linear increase with the sediment grain size in the point bars ( Figure 6). Vertically, the region where the permeability begins to increase rapidly corresponds approximately to the middle section of the sandstone. The vertical permeability distribution function is introduced as an auxiliary condition of the permeability simulation in fluvial channels. In contrast, the permeability distribution function in the lateral accretion direction is employed in the permeability simulation in point bars (Figure 14(d)). Moreover, different exponential functions are applied in the different point bars. Finally, the minimum truncated value of the permeability of the argillaceous intercalations was set at 0.1 mD in the reservoir architectural model (Figure 14(c)).

Model Verification
The goal was to evaluate whether this integrated approach to reservoir characterization, which includes constructing the depositional model and multiscale 3D geomodeling of point 11 Geofluids bars and their internal heterogeneities, is practical. The connected-well section parallel to the orientation of the lateral accretion layers was selected for visual inspection in the dense well network area ( Figure 15). The proportion of the sedimentary facies was consistent in the simulation realization, the well data, and the TIs (Figure 16(a)).
Dynamic production data can also be used for model verification because the distribution of reservoir physical properties affects historical production. A reservoir numerical simulation is based on property models, in which the liquid product in the model area is fitted with the actual data. As a result, the water cut rate of the MPS realization is sim-ilar to the actual data, which partially reflects the reliability of the simulation (Figure 16(b)).
Furthermore, two wells were selected to verify that the lateral accretion layers affect the horizontal permeability. Well Zh8-13 and Well Zh9-13-3 are a pair of wells in the same point bar with a distance of 100 m (Figure 14(b)). After a production period, Well Zh8-13 stopped production and transformed into a water injection well due to a high water-cut rate. However, due to the effect of the lateral accretion layer barrier, the injected water barely affected the displacement of oil in the upper section sand body of Well Zh9-13-3. This may be due to a high permeability zone   13 Geofluids and a dominant seepage channel formed at the bottom of the reservoir (Figure 15). This result also verified the accuracy of the reservoir model. 14 Geofluids are composed of clast-supported conglomerate, massive structureless sandstone, cross-bedded sandstone, horizontally bedded sandstone, ripple laminated siltstone, horizontally laminated siltstone, and massive mudstone, with an associated fining-upward grain size profile. As an impermeable barrier screen, the lateral accretion layers of mudstone mainly extend from the top to 1/2 or 2/3 the thickness of the point bar. These lateral accretion layer data identified by electrofacies data can be used to improve the reliability of the multiscale facies model containing lateral accretion layers

Conclusions
(3) The permeability of the lateral accretion elements exhibits exponential variation with depth. The vertical distribution function of permeability is treated as an auxiliary variable for collaborative simulation of the channel permeability model. In comparison, the permeability model of point bars is simulated collaboratively by the permeability distribution function in the direction of the lateral accretion elements. Finally, model verification by three methods illustrates the accuracy and practicability of the model

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request. 15 Geofluids