The reliability of forecasts, fracture design, and recovery enhancement strategies in tight oil reservoirs is significantly compromised by the substantial uncertainties associated with fracture characterization. This article introduces an integrated simulation workflow for modeling microseismic fracture networks in tight oil reservoirs, incorporating automatic history matching, as illustrated through a field case study from block Y2 in the Ordos Basin, China. The model stochastically generates the geometry of complex fracture networks (CFNs), including parameters such as length, aperture, inclination and azimuth angles, and spatial positioning, constrained by data from hydraulic fracturing, core analyses, and microseismic monitoring. It employs a stochastic parameterization model to produce an ensemble of initial CFN property realizations and utilizes an advanced Green function-based hierarchical fracture model to accurately depict CFN morphology. The model is further refined, and its uncertainty in fracture characterization is minimized through calibration with an innovative Ensemble Kalman Filter-based assisted history-matching algorithm. Evidence suggests that this comprehensive approach effectively leverages all available geological data, substantially reduces uncertainties in the production process, and aids in identifying the optimal development strategy.

Unconventional reservoirs are characterized by ultra-low permeability, high displacement resistance, and low productivity, primarily due to the presence of nanopores and fine throats within the tight matrix. Extensive experience in successfully extracting hydrocarbons from these reservoirs has demonstrated that creating complex fracture networks (CFNs) through multistage hydraulic fracturing represents one of the most effective strategies for boosting oil production [1-3]. However, accurately modeling CFNs presents significant challenges due to substantial uncertainties, leading to increased difficulties in the precise numerical simulation of unconventional reservoirs [4]. Additionally, the disconnect between geological and petroleum engineering disciplines impedes reservoir engineers from fully leveraging geological data in production planning, constituting another significant factor contributing to these challenges.

Obtaining fracture morphology through conventional characterization methods is challenging. To enhance the assessment of fracture occurrence, propagation, and the effects of refracturing, microseismic monitoring has been employed to provide essential geological information, as illustrated in Figure 1. While microseismic monitoring data, collected during hydraulic fracturing, enables the analysis of fracture spatial distribution, accurately characterizing CFNs from microseismic events remains a technically demanding task. Currently, methods for reconstructing fracture networks based on microseismic interpretation are categorized into two approaches. The first relies on hydraulic fracture propagation models (such as PKN, KGD, and P3D), employing either the displacement discontinuity method [5-7] or finite element method [8-10] in conjunction with microseismic data for dynamic simulation of fracture propagation. However, this method overly depends on the modeling process, resulting in underutilization of microseismic data. The second approach uses a linear relationship to approximate the trend of microseismic points for constructing primary fracture networks [11, 12]. Although quick and straightforward, this method is susceptible to various subjective biases from reservoir engineers and often does not correlate well with actual well production data.

Characterizing fractured reservoirs numerically is notably challenging. The prevalent approaches include the classic multicontinuum model, exemplified by DPDP or MINC, and explicit discretization methods, such as the Discrete Fracture Model (DFM) [13-15] and the Embedded Discrete Fracture Model (EDFM) [16-18]. Multicontinuum models often oversimplify reality, whereas DFM incurs high computational costs due to fine mesh generation. EDFM, on the other hand, combines the benefits of fine mesh refinement with mesh simplicity and flexibility, significantly reducing preprocessing computational demands. Consequently, EDFM has garnered considerable attention and advanced significantly within the academic sphere. Jiang et al. [19, 20] extended EDFM to nonorthogonal unstructured grids with a mimetic finite difference approach and developed a multicomponent model for analyzing gas injection in shale reservoirs. Similarly, Yu et al. [21] utilized EDFM to study the impact of shale gas transport on the production performance of multistage fractured horizontal wells. Despite its successful applications across various domains, EDFM still faces unresolved fundamental challenges. Rao et al. [22, 23] integrated the concept of Green elements into EDFM, establishing a mimetic Green element method (GEM) that approximates the boundary integral in Green elements for normal flux calculation in finite differences. This innovation has led to the creation of next-generation simulators that are highly accurate and feature enhanced functionality for 3D EDFM coupling. Moreover, simulation studies reveal that oil reservoir production performance is highly sensitive to fracture parameters like geometry, connectivity, and conductivity. Nonetheless, initial fracture modeling is fraught with uncertainties, stemming from the vagueness of microseismic data and descriptive methods. The complexity and burdensome nature of history-matching processes render manual operations inefficient and prone to subjective bias. Hence, there is a pressing need for the development of computer-aided history-matching technologies to mitigate human error and enhance simulation accuracy.

Automatic history matching, or parameter inversion in oil reservoirs, is a technique aimed at optimizing the objective function to its minimum value. This process involves various objective functions, with differences among parameter inversion methods stemming from the selection of optimization algorithms to define their specific forms. The methods to achieve this optimization are categorized into gradient-based and nongradient approaches. Gradient-based methods, noted for their simplicity and effective convergence, have seen widespread application. For instance, Bissell et al. [24] applied a direct-solving gradient method to assess sensitivity coefficients in history matching, while Chen [25] and Chavent [26] leveraged the adjoint equation of gradient methods for optimal control problems. However, these methods face challenges in large-scale applications due to the extensive number of model parameters involved. Nongradient methods, including evolutionary algorithms, artificial neural networks (ANNs), and Bayesian theory, offer alternative strategies. Evolutionary algorithms, distinguished by their stability and broad applicability, have been successfully applied by Perez et al. [27] to enhance simulation efficiency in fractured reservoirs, integrating with EDFM and accounting for water channeling effects, yielding more reliable outcomes than traditional approaches, despite issues with slow convergence [28]. ANNs, capable of adapting to diverse scenarios and identifying optimal solutions through learning, often require integration with other techniques, adding complexity. The Ensemble Kalman Filter (EnKF), based on Bayesian theory, excels in updating parameter models with new data, maintaining the predictive capacity of historical data, and significantly streamlining history-matching processes [29, 30]. Its application is particularly beneficial in fractured reservoir models, where it can enhance the efficiency and accuracy of numerical simulations by adjusting model parameters, such as fracture geometry. Liu and Dai’s combination of EnKF with DFM for estimating fracture distribution and Wu et al.’s use of EnKF with EDFM to minimize uncertainties in fracture and reservoir parameters of multistage fractured horizontal wells exemplify its superiority in fractured reservoir simulations [31, 32]. Nonetheless, the approach generally focuses on single-phase fluids and struggles with accurately interpreting microseismic fracture networks, highlighting a gap in practical applications and case studies in tight oil reservoirs.

An efficient workflow integrating techniques such as microseismic fracture modeling, numerical simulation, and calibration through assisted history matching has been developed for a case study in the Ordos Basin, China. Initially, an effective forward simulation method was utilized to preliminarily characterize fracture networks based on microseismic event points. Subsequently, an enhanced Green function-based EDFM was employed for accurately modeling the morphology of hierarchical fractures. Finally, the EnKF was applied for computer-aided history matching, identifying sensitive parameters that influence various performance indicators and determining how to adjust these variables for precise production forecasting. This approach effectively addresses the challenge of accurately calibrating and characterizing microseismic fracture networks, a task that conventional methods fail to achieve. The field case simulation results affirm the workflow’s feasibility and efficiency.

2.1 Geological Context of Study Area in Ordos Basin

The Ordos Basin, situated within the transitional zone between the eastern and western tectonic units on the North China Platform, is depicted in Figure 2. This basin is characterized by its multicycle cratonic nature, exhibiting stable subsidence during the Paleozoic era, a west-to-east migration of Mesozoic depressions, and significant twisting and faulting throughout the Cenozoic period [33]. The area of interest, Block Y2 in the HQ Oilfield, is located in Wuqi County, Gansu Province, within the southwestern part of the Yishan slope in the Ordos Basin. This region features a typical loess plateau topography and a gentle monocline structure with formation dip angles below 1°. Despite several small-scale nose uplifts, the area lacks significant fault development. The primary oil-bearing layer, the Chang-6 of the Triassic period, is classified within lithological reservoirs and is primarily influenced by a marine sedimentary environment. This layer has an average thickness of 18.5 m and exhibits pronounced vertical heterogeneity. Matrix characteristics are predominantly marked by micropores and fine throats, typifying it as a tight reservoir. Challenges such as high seepage resistance, difficult displacement, and low recovery efficiency during development have necessitated the adoption of well-factory models employing horizontal wells and extensive hydraulic fracturing for the efficient exploitation of tight oil reservoirs.

2.2 Fracture Identification and Analysis

The characteristics of natural fractures in Block Y2 are extensively characterized through conventional logging and core observation, as illustrated in Figure 3. These fractures are prevalent and abundant, with fracture densities exceeding 0.6 fractures/m in 50% of the analyzed wells, over 1 fracture/m in 25% of the wells, and an overall average density surpassing 0.76 fractures/m. The fractures generally exhibit small apertures, with approximately 90% of the wells showing apertures narrower than 0.13 mm, and around 60% displaying apertures under 5 mm. Furthermore, the penetration depth of these fractures is limited, with about 75% not exceeding 20 cm, indicating that they typically have narrow apertures, close spacing, shallow penetration, and frequent fillings. Examination of the dip angles of core fractures reveals that the block predominantly features high-angle fractures, with dip angles averaging 87°. Comprehensive imaging logging and outcrop analysis have determined the orientation of natural fractures in Block Y2, as depicted in Figure 4, to be multidirectional, predominantly extending from NE 35° to 75°, with a general NE to SW orientation. Core analysis of stress parameters indicates a relatively minor horizontal stress differential, ranging between 5 and 6 MPa, suggesting the feasibility of developing CFNs.

To enhance the optimization of hydraulic parameters within fracture networks, microseismic monitoring was conducted on two multistage fractured horizontal wells, CP-1 and CP-2, as depicted in Figure 5. This monitoring enabled the acquisition of spatial distribution data for fractures during hydraulic fracturing processes. For CP-1, 22 fracturing stages were implemented, achieving a total reservoir reconstruction volume of approximately 4909.9 × 104 m3. The microseismic results indicated a notable disparity in fracture lengths between the wings, with the west wing fractures extending significantly beyond 200 m—markedly longer than the east wing’s fractures, which predominantly fell within 150 m. This variation primarily arises from reservoir heterogeneity. Similarly, CP-2 underwent fracturing in 22 stages, reaching a total reconstruction volume of about 7348.9 × 104 m3. Notably, the initial sections (first to fourth) displayed longer fractures on the east wing compared to the west, whereas fractures from the fifth to the twenty-second section showed more uniform lengths between wings. Figure 6 outlines the workflow for fracture identification and analysis within reservoirs, incorporating both natural and hydraulic fractures, though the latter is not detailed here.

Determining key parameters for fractured horizontal wells presents substantial challenges due to the complex nature of effective fracture length, attributes of natural fractures, and hydraulic fracture conductivity among others. Addressing the difficulty in acquiring these uncertain parameters is crucial. Utilizing existing datasets for accurate fracture characterization emerges as a viable strategy to mitigate overall parameter uncertainty. This study leverages a diverse range of data—including core, microseismic, hydraulic fracturing, and production performance data—as constraints for detailed fracture analysis. It outlines a comprehensive procedure for parameter estimation, incorporating techniques such as the microseismic fracture forward modeling, an enhanced GEM-based EDFM for numerical simulations, and the EnKF approach for assisted history matching.

3.1 Microseismic Fracture Forward Modeling

The primary objective in analyzing microseismic data is to infer critical parameters of natural fractures. Mayerhofer et al. [34] suggested that natural fractures within reservoirs reopen due to hydraulic forces, with the re-tension of fractures releasing energy as seismic waves, thereby generating microseismic events. Consequently, the coordinates of these microseismic events serve as proxies for the locations of the re-tensioned natural fractures. To model natural fractures accurately, three essential parameters are identified: the location, length, and orientation of the fractures [35]. The location is determined by microseismic event coordinates, assuming each event corresponds to a single natural fracture. Furthermore, core and image log analyses provide the basis for estimating the orientation, length, and abundance of natural fractures within the reservoir, facilitating the statistical distribution of these characteristics. With these relational parameters established, the length and direction of each fracture are generated randomly based on the observed fracture distribution. This framework enables the simulation of both the spread of hydraulic fractures through the matrix and their interaction with existing natural fractures. Olson [36] demonstrated through experimental physics that in the absence of natural fractures, hydraulic fractures propagate in the direction of the maximum horizontal stress. The behavior of hydraulic fractures upon intersecting with natural fractures significantly depends on the angle of intersection and the resulting stress-induced displacement.

The aforementioned research has provided significant insights. In prior studies [37], we introduced an enhanced EDFM to simulate the development of hydraulic fractures, achieved by refining the preprocessing algorithm. This methodology facilitates fracture growth simulation by appending new fracture elements to existing ones upon meeting the failure criterion. To adequately address the effects of fracture interactions across various stages, we treat hydraulic fractures from earlier stages as preexisting fractures prior to the propagation in subsequent stages. Thus, hydraulic fractures evolve progressively from the wellbore’s toe, interacting with existing fracture networks. Our approach addresses the issue of unlimited fracture propagation seen in traditional models by positing that fracturing fluids swiftly reactivate natural fractures and follow their path until intersecting with another natural fracture, halting the hydraulic fracture’s extension when it attains the predetermined length for that stage. Furthermore, we hypothesize that the total length of fractures created during each fracturing phase is directly proportional to the volume of fracturing fluids utilized, thereby establishing a quantitative relationship between fluid volume and fracture extension, as shown in equation. (1).

Xf,i=βfVf,i,
(1)

where Xf,i represents the length of hydraulic fractures generated for the ith stage; βf represents a constant coefficient; and Vf,i represents the volume of fracturing fluids in the ith stage.

A crucial aspect involves identifying propped fractures and quantifying the length of proppant-filled fractures. Through a combination of experimental research and numerical simulations, Tong et al. [38] demonstrated that when hydraulic fractures intersect with natural fractures, only a fraction of the proppant manages to enter the natural fractures, leaving the majority within the hydraulic fractures. Consequently, our research posits that proppant is predominantly present in primary fractures, with interconnected induced fractures lacking proppant support. Additionally, we hypothesize that the total length of proppant-filled fractures bears a direct proportionality to the overall quantity of proppant used, as illustrated in equation. (2).

Xp,i=βpmp,i,
(2)

where Xp,i represents the total length of level i proppant-filled fracture; βp represents a constant coefficient; and mp,i represents the total amount of proppant used for fracturing stage i.

This study illustrates the visualization of fracture propagation simulation using microseismic data from a specific section of the CP-1 fractured horizontal well. Figure 7(a) presents the microseismic event map along with pertinent details. The analysis was conducted on a grid with a 100 m resolution. The total injected liquid and sand volumes were 1857.68 and 120 m3, respectively. This data enables the identification of potential fracture initiation points in the reservoir rock, facilitating the creation of a coordinate system to simulate natural fractures at microseismic event locations. The distribution of fracture lengths follows a normal distribution, with an average length of 40 m and a standard deviation of 20 m. Assuming the known fracturing perforation point as the origin, hydraulic fractures extend bilaterally, perpendicular to the wellbore, under the influence of the maximum principal stress direction. This process leads to encounters with natural reservoir fractures, either bypassing or integrating with them to form CFNs, as depicted in Figure 7(b). Given the computational demands of simulating free-propagating fracture morphology, a simplification to represent fractures as straight lines in the numerical model was employed, where red signifies propped hydraulic fractures and blue indicates natural fractures without proppant. This simulation exemplifies the feasibility of the forward-modeling approach detailed in this study.

Upon determining the geometry of fractures, it becomes essential to define their properties, specifically focusing on aperture and conductivity. The distinction between propped and unpropped fractures necessitates assigning their properties separately due to the significant variance in their characteristics. Barton et al. [39] introduced a joint-closure relationship, represented as equation (3), to estimate the aperture of closed segments in both propped and unpropped fractures. A similar formula [40], with alternate constants for calculating hydraulic aperture, is described as equation (4). This study employs these nonlinear equations to delineate the properties of fractures individually. Furthermore, adhering to the classical cubic law proposed by Witherspoon et al. (1979) [41], fracture conductivity is correlated with aperture, as outlined in equation (5). These equations serve as the basis for generating an initial set of fracture parameters. Subsequent refinement of fracture properties is achieved through history matching.

w=w01+9σcσc,wref,
(3)
wF,prop=w01+9σcσc,wref+Dw,efftanφwdil1+9σcσc,wref,
(4)
kFw=w312,
(5)

where w represents the fracture aperture; w0 is the contact aperture when the contact stress is equal to 0; σc represents the contact stress between the fracture surface and the rough surface of proppant packs; σc,Wref represents the contact reference stress, which denotes the effective stress at which the fracture aperture is reduced by 90%; wF,prop is the hydraulic aperture; φWdil represents fracture dilation angle; and kF represents the fracture permeability.

3.2 Review of Previous Improved GEM-Based EDFM

In our prior research [42, 43], we developed a Hierarchical Fracture Model (HFM) utilizing an enhanced GEM to delineate fracture characteristics across various scales and to accurately simulate fluid flow in porous media. The refinement lies in employing a modified EDFM to distinctly represent macro-scale hydraulic fractures, while small-scale natural fractures within the Stimulated Reservoir Volume (SRV) are approximated using a dual-medium model. Thus, the advanced GEM serves as an effective augmentation of the boundary element method, adept at addressing nonlinear and heterogeneous challenges in fractured reservoirs. This novel approach surpasses traditional Finite Difference Method (FDM)-based EDFM simulators in precision, particularly in capturing flux transient behaviors during initial flow periods. Figure 8 illustrates the methodology for managing flux in EDFM within a dual-porosity system. This fracture modeling strategy will be further applied in the current study.

3.3 The EnKF Approach for Automated History Matching

Previously discussed, EnKF serves as an auxiliary method for history matching, particularly effective for numerical simulations in fractured reservoirs. Within the EnKF framework, an ensemble of initial models, derived from prior information, is utilized for predictive simulations. These models are subsequently updated upon the acquisition of new data. The relationship between the updated data and the model’s state vector is established through statistical analysis of the ensemble. This technique is notably advantageous for addressing complex challenges in reservoir simulation [44]. To delineate the algorithm, we first define the state vector as follows:

y˜k={muk}RNm+Nu,
(6)

where y˜k is the state vector; m is the static variable of the model, such as porosity and permeability; u is the dynamic variable, such as the pressure and saturation in reservoirs; the subscript k is the time step; R is the real vector space; and Nm and Nu are the dimensions of vectors m and u, respectively.

In the sequential method, it is necessary to discuss conditional probability p(y˜kDobs,k). According to Bayesian theory, then:

p(y˜kDobs,k)=p(m,ukDobs,k)p(dobs,km,uk)p(m,ukDobs,k1)=p(dobs,ky˜k)p(y˜kDobs,k1),
(7)
Dobs,k={dobs,1,dobs,2,dobs,k},
(8)

where dobs,k represents the observation value of production data and Dobs,k represents the observation value vector of production data.

If it is assumed that both the overall error and the state vector y˜k obey the joint Gaussian distribution, then

p(y˜kDobs,k)=cexp{12[g(y˜k)dobs,k]TCDk1[g(y˜k)dobs,k](y˜ky˜kp)TCy˜kp1(y˜ky˜kp)},
(9)

where CD is the covariance matrix of measurement errors of production data; g(y˜k) is the predicted value of production data, Cy˜kp is the sensitivity gradient matrix of dynamic production measurement; the superscript P represents the predicted value; and c is a constant.

The objective function is

J(y˜k)=12[g(y˜k)dobs,k]TCDk1[g(y˜k)dobs,k]+12(y˜ky˜kp)TCy˜kp1(y˜ky˜kp)
(10)

where J(y˜k) is the objective function.

Using the full-step Gauss–Newton iteration, and then:

y˜k,ja=y˜k,jp+CY˜kpGk,jT(CDk+Gk,jCY˜kpGk,jT)1[dobs,k,jg(y˜k,jp)],(j=1,2,,Ne),
(11)

where y˜k,j is the state vector corresponding to the jth reservoir model; dobs,k,j is the observed value after disturbance; Gk,j is the sensitivity coefficient matrix; Ne is the number of reservoir models; and the superscript a indicates update.

The right end of the equation (11) can be expressed as [45]:

CY˜kpGk,jT(CDk+Gk,jCY˜kpGk,jT)1CY˜kpG¯k,jT(CDk+G¯k,jCY˜kpG¯k,jT)1CY˜kp,Dk(CDk+CDk,Dk)1,
(12)

where G¯k represents sensitivity parameters; G¯k,j represents the sensitivity matrix; CY˜kp,Dk represents the interaction covariance matrix of model state vector predicted values Y˜kp and dobs,k; and CDk,Dk is the autocovariance matrix of predicted values of production data.

Therefore, it can be expressed as:

y˜k,ja=y˜k,jp+CY˜kp,Dk(CDk+CDk,Dk)1[dobs,k,jg(y˜k,jp)],
(13)

To achieve sequential algorithms, the state vector of the model needs to include the variable u. In EnKF, model state vectors include the parameter vector m, the state variable u, and the observation data d, as:

y˜k={mukdk}RNm+Nu+Nd,
(14)
dk=Hkyk,
(15)
Hk=[OI],
(16)

where O represents Nd×(Nm+Nu)-dimensional matrix, all elements in matrix are 0; Nd is the spatial dimension; Hk is the observation matrix; yk is the observation state vector; and I represents Nd×Nd-dimensional identity matrix.

After introducing the observation matrix, we can further obtain that

CYkp,Dk=CYkpHkT,
(17)
CDk,Dk=HkCYkpHkT,
(18)

where CYkp is the covariance matrix of state vectors and CYkp,Dk is the interactive covariance matrix of predicted values of state vectors of dynamic production gradient matrix model and production data predicted values.

Using equations (13), (17), and (18), the updated expression of EnKF can be obtained:

yk,ja=yk,jp+CYkpHkT(CDk+HkCYkpHkT)1(dobs,k,jHkyk,jp),
(19)

where yk,j represents the state vector corresponding to the jth reservoir model.

EnKF includes two parts: prediction and update. During the predicting process, all data advances in time independently of each other, as:

yk,jp=F(yk1,ja),
(20)

where F is the prediction operator and also the reservoir simulator.

The Kalman kernel matrix Kk is defined as equation (21), then the update formula can be transformed into equation. (22).

Kk=CYkpHkT(HkCYkpHkT+CDk)1,
(21)
yk,ja=yk,jp+Kk(dobs,k,jHkyk,jp).
(22)

The EnKF process involves the following steps: Initially, users must employ statistical techniques to create datasets based on available empirical data or information about the physical field. The quantity of generated data typically should be double that of the variables in question. Once this data generation phase is complete, parameters such as permeability are updated through EnKF. The timing for these updates is set based on user definitions and the locations of the generated data. EnKF facilitates not only the modification of parameters but also the adjustment of state vectors, including saturation and pressure levels. Prior to proceeding with subsequent computational steps using these new parameters, it is imperative to refresh the currently calculated saturation and pressure fields. This action effectively recalibrates the initial conditions for the forthcoming calculations, as depicted in Figure 9. This iterative process is continued until all updates have been applied, culminating in a parameter field that accurately reflects the underground reservoir’s physical characteristics and fluid dynamics to the fullest extent.

This study focuses on horizontal wells CP-1 and CP-2 within Block Y2 of the HQ Oilfield, examining their refracturing processes. The tight reservoirs exhibit an average matrix porosity of 0.12 and an average matrix permeability of 0.34 millidarcies (mD). Both wells underwent refracturing, with CP-1 being segmented into 22 stages, where each stage involved the injection of 1897 m³ of fracturing fluids and 134 m³ of proppant. Similarly, CP-2 was fractured through 23 stages, each receiving an average of 1683 m³ of liquid fluids and 134 m³ of proppant for injection. Microseismic monitoring spanned sections 1–22 for CP-1 and sections 2–23 for CP-2, facilitating the reconstruction, simulation, and refinement of active fracture networks employing the methodologies discussed earlier.

4.1 Reconstruction of Microseismic Fracture Networks

Identifying viable microseismic data points within the study area and discarding unreliable data involves two key steps: (1) Energy Filtering: isolating points where microseismic energy falls within an interpretable and reasonable range; (2) Spatial Filtering: excluding distant microseismic points based on their spatial location during the sequential analysis of fracture network reconstruction. The modeling of natural fractures necessitates the definition of primary parameters including the location, length, and orientation of fractures. Microseismic events determine the locations of natural fractures, with the assumption that each event generates a single natural fracture. The height of these natural fractures aligns with the reservoir’s thickness and is positioned at the reservoir’s midpoint. Core and image log analyses provide data on the azimuth and inclination angles, length, and aperture, enabling the establishment of a frequency distribution for the attributes of each fracture, as listed in Table 1. Subsequently, natural fractures are randomly generated and modeled according to this distribution.

The characteristics of hydraulic fractures are determined by pumping properties, with the pumping capacities of CP-1 and CP-2 wells detailed in Tables 2 and 3, respectively. The microseismic fracture network is modeled by simulating the interaction and expansion of hydraulic and natural fractures, as depicted in Figure 10. Among them, Figure 10(a) depicts the location of microseismic events as black circles, with lines within these circles indicating the orientation of natural fractures. Figure 10(b) presents the varied scales of fractures: hydraulic fractures are marked by red lines, induced fractures by blue, and natural fractures by green. This imagery reveals that volume fracturing in the tight reservoirs of Block Y2 generates intricate fracture networks characterized by bilateral extension and localized branching. Such formations are significantly influenced by factors including pumping parameters, in situ stress levels, and the prevalence of natural fractures. Although microseismic networks signal the presence of subterranean fractures, not all exhibit effective conductivity, necessitating further analysis for refinement. By extracting valid fracture data from Figure 10—referring to actively functioning fracture networks post-elimination of spurious microseismic modeling data, as demonstrated in Figure 11—the effective lengths of hydraulic fractures are highlighted in red. Post-refracturing, the semi-length of these fractures ranges from 65.3 to 260.3 m, averaging at 141.7 m. This indicates that refracturing notably enhances the structural complexity of fracture networks.

4.2 History Matching and Numerical Simulation

This study utilized an assisted history-matching method anchored in real production data to refine the delineated effective fracture networks. History matching represents a complex process of multiparameter inversion, necessitating tight constraints on the variability of each parameter to prevent distortions, particularly for those parameters that exhibit high uncertainty and significantly influence predictive outcomes. Identifying and constraining the variability of these critical parameters ensures the numerical simulation results are both credible and robust. Among these parameters, reservoir properties such as porosity, permeability, and the oil–water interface are inferred from a detailed analysis of log data. Meanwhile, pressure–volume–temperature (PVT) data and relative permeability curves for reservoir fluids are derived from laboratory experiments, underscoring the importance of evaluating the reliability and sensitivity of these parameters. For instance, assessing the validity of the water phase permeability curve or the reliability of initial water saturation becomes imperative when actual water production exceeds model predictions. In fractured reservoirs, the complexity and heterogeneity of fracture networks—characterized by their spatial distribution, porosity and permeability, and fracture aperture—demand cautious evaluation and potential adjustment during history matching. Geomechanical properties and rock physical parameters also present considerable uncertainties and risks. The multiscale fractures depicted in Figure 11 are effectively modeled using a HFM with dual-medium grids, as illustrated in Figure 12. A comprehensive assessment of key sensitivity parameters through the EnKF automated history-matching process is detailed in Table 4.

Production data from wells CP-1 and CP-2, spanning the initial fracturing phase from December 2012 to May 2018 and the refracturing phase from June 2018 to August 2021, were utilized for assisted history matching. The refined outcomes, following several iterations of updates, are depicted in Figures 13 and 14. These outcomes reveal that the modeled daily oil production closely aligns with the actual data, maintaining a margin of error within 5%. A comparison of the effective fracture half-lengths derived from history matching and the forward-modeled microseismic fracture network outcomes, as presented in Tables 5 and 6, underscores disparities between the calculated effective fracture half-lengths and the microseismic observations. These discrepancies range from 2% to 6% across various parameters, highlighting the importance of critically evaluating forward-modeling results of identified effective fracture networks. Leveraging the EnKF for automated history matching notably diminishes the inaccuracies associated with subjective interpretations, thereby enhancing the precision of fracture network parameters. By conducting automated history matching with production data, the reconstructed effective fracture networks have been affirmed to satisfy simulation accuracy requirements, thus offering a dependable framework for forecasting future production.

Building on the foundation of reliable parameters for CFNs, a numerical model integrating dual-porosity media with embedded discrete fractures, as illustrated in Figure 12, is employed for development index prediction. The enhanced GEM-based EDFM is utilized to represent propped hydraulic fractures and induced fractures, while the dual-porosity media models the SRV area filled with activated natural fractures. This model adeptly leverages corrected reservoir parameters, fluid properties, and fracture network details to forecast production in CP-1 and CP-2 wells under two distinct development scenarios: (1) advanced water injection followed by refracturing, and (2) development without energy supplementation measures. Figure 15 displays the predicted pressure distribution maps across various time frames. Figure 16 depicts the daily oil production trend throughout the entire cycle, encompassing the initial fracturing, advanced water injection, well shut-in, refracturing, and production forecast phases. The production forecasts are generated through numerical simulations based on the automated history matching of production data for wells CP-1 and CP-2. The results indicate that following advanced water injection and refracturing, the anticipated subsequent production for both wells markedly exceeds that of scenarios lacking energy supplementation. Thus, employing water injection to augment formation energy prior to refracturing proves to be an effective strategy for enhancing productivity in tight oil reservoirs.

This article introduces an integrated workflow combining forward modeling and inversion correction for a field case, offering several benefits. Primarily, forward modeling of CFNs utilizing microseismic data accurately depicts the distribution and expansion of hydraulic fractures. Advanced fracture modeling employing the GEM alongside structured grids enables precise characterization and simulation of fluid dynamics in tight oil reservoirs. The oil production estimates, derived through EnKF-assisted history matching, align closely with actual field data. By merging reservoir simulation, microseismic fracture forward modeling, detailed fracture network characterization, and EnKF-assisted matching, this approach effectively harmonizes historical production data with clear delineation of the CFNs system and minimized uncertainty. Consequently, it facilitates the precise prediction of oil production in tight reservoirs.

No additional data are available.

The authors declare that they have no conflicts of interest.

This research received funding from the National Natural Science Foundation of China under grant number 52174038, and the China Petroleum Science and Technology Program, specifically the Major Program ZLZX2020-02-04.

Thank you to the professors and students of the Fluid Flow Mechanics and Reservoir Numerical Simulation Team at the College of Petroleum Engineering, China University of Petroleum (Beijing) for their help and support. In particular, I would like to express my gratitude to Dr. Shi Junjie and Dr. Yang Chenxu for their academic guidance.

Exclusive Licensee GeoScienceWorld. Distributed under a Creative Commons Attribution License (CC BY 4.0).