A method for establishing a compact reservoir matrix fracture enzyme diffusion model

Abstract

The application provides a method for establishing a compact reservoir matrix fracture enzyme diffusion model, and belongs to the technical field of compact reservoir matrixes.The method comprises the following steps: collecting a compact reservoir core sample; constructing a multi-scale effective transfer coefficient matrix of hierarchical transmission by using a multi-scale homogenization theory; combining a discrete fracture network random generation technology with an integrated Kalman filter data assimilation method to constrain fracture parameter posterior probability distribution and output an equivalent permeability tensor and a matrix fracture mass exchange function; and dynamically adjusting a physical information loss weight coefficient by using a concentration field self-adaptive regulation function to drive iteration until a global error tolerance is satisfied, and finally outputting a matrix fracture enzyme diffusion model parameter set, so as to solve the technical problem that, in a compact reservoir matrix fracture double-pore medium, a multi-scale mass transfer coupling calculation of an enzyme is insufficient in physical self-consistency in time sequence prediction of a cross-scale concentration field due to fracture topological uncertainty.

Description

Technical Field

[0001] This invention belongs to the field of tight reservoir analysis technology, and more specifically, relates to a method for establishing a biological enzyme diffusion model of matrix fractures in tight reservoirs. Background Technology

[0002] Tight reservoir bio-enzyme displacement technology, which involves injecting ESS-100 and ELS bio-enzymes into the reservoir, utilizes enzymatic reactions to improve pore permeability and wettability, and is an important production enhancement method in unconventional oil and gas development. In traditional modeling practices, researchers typically use the Darcy flow equation at a single scale to describe the migration process of bio-enzymes in the reservoir, characterize the fracture network using a deterministic fracture geometry model, and numerically solve the concentration field distribution using finite difference or finite element methods. This approach has achieved certain engineering applicability in homogeneous reservoirs or low fracture density scenarios.

[0003] However, in tight reservoirs, three mass transfer mechanisms coexist: Knudsen diffusion in nanopores, Fick diffusion in micropores, and convection mass transfer through fractures. Single-scale models cannot coordinate the governing equations of these three mechanisms within a unified mathematical framework, leading to systematic errors in the mass transfer coefficient across scales. Simultaneously, the intrinsic randomness of subsurface fracture geometry and topology causes ill-conditioned Jacobian matrices in high-dimensional parameter inversion using deterministic fracture models, resulting in significant deviations of the equivalent permeability tensor estimation from the true reservoir state. Furthermore, traditional recurrent neural networks suffer from the vanishing gradient problem in long-term concentration field prediction and lack physical constraints, making it difficult to maintain physical consistency of prediction results under sparse labeling conditions. In other words, existing technologies suffer from technical problems such as ill-conditioned equations due to fracture topological uncertainties and insufficient physical consistency in cross-scale concentration field time-series predictions in the dual-porosity medium of tight reservoir matrix and fractures. Summary of the Invention

[0004] In view of this, the present invention provides a method for establishing a biological enzyme diffusion model in tight reservoir matrix fractures, which can solve the technical problems in the prior art where the multi-scale mass transfer coupling calculation of biological enzymes in tight reservoir matrix fracture dual-pore media is ill-conditioned due to fracture topological uncertainty, and the physical self-consistency of cross-scale concentration field time series prediction is insufficient.

[0005] This invention is implemented as follows: This invention provides tight reservoir analysis, including the following steps:

[0006] Core samples from tight reservoirs were collected, and pore-scale distribution parameters from nanopores to millimeter-scale fractures were obtained through nitrogen adsorption experiments and mercury intrusion porosimetry experiments. Based on the multi-scale homogenization theory, the Knudsen diffusion coefficient of nanopores, the Fick diffusion coefficient of micropores, and the fracture convection mass transfer coefficient were asymptotically expanded and uploaded to the Darcy scale to construct a multi-scale effective transport coefficient matrix of hierarchical transfer.

[0007] The discrete fracture network random generation technique is used to construct integrated samples of fracture geometric topology for pore-scale distribution parameters. Dynamic pressure observation data is injected into the Monte Carlo forward propagation process through integrated Kalman filter data assimilation method to constrain the posterior probability distribution of fracture aperture, fracture orientation and fracture connectivity, and output the equivalent permeability tensor and matrix fracture mass exchange function.

[0008] Using the equivalent permeability tensor and matrix fracture mass exchange function as input, the matrix fracture coupled diffusion reaction equation is decomposed into diffusion substeps and enzyme-catalyzed reaction substeps. An adaptive multigrid enzyme kinetics operator splitting acceleration algorithm is used to solve the problem iteratively on a multi-grid. After each complete operator splitting cycle, the coupling residual is calculated and the number of local grid layers and the number of smoothing iterations are dynamically adjusted. The numerical solution of the grid node concentration field at each time step is output.

[0009] The multi-scale effective transport coefficient matrix, the equivalent permeability tensor and the matrix fracture mass exchange function, and the numerical solution of the concentration field of the grid nodes are input into the enzyme concentration spatiotemporal evolution prediction model. The enzyme concentration spatiotemporal evolution prediction model outputs the three-dimensional concentration field time series prediction results of ESS-100 bioenzyme and ELS bioenzyme in the matrix fracture dual-porosity medium.

[0010] Using the coupling residual, the posterior standard deviation of the crack parameters, and the prediction residual of the three-dimensional concentration field time-series prediction results as inputs, the value of the concentration field adaptive control function is calculated. The physical information loss weight coefficient of the enzyme concentration spatiotemporal evolution prediction model is dynamically adjusted according to the interval to which the value of the concentration field adaptive control function belongs, and the adaptive multigrid enzyme kinetic operator splitting acceleration algorithm is re-driven to iterate until the global error tolerance is met.

[0011] Based on the convergent three-dimensional concentration field time-series prediction results that meet the global error tolerance, and combined with the statistical results of the posterior probability distribution of crack parameters, the spatial distribution of each component of ESS-100 bioenzyme and ELS bioenzyme in the matrix and cracks is evaluated, and the final matrix crack bioenzyme diffusion model parameter set is output.

[0012] Specifically, the multi-scale homogenization theory uses an asymptotic expansion method as a mathematical tool to eliminate the microscopic heterogeneity in the pore-scale control equations through periodic cell integration, obtaining macroscopic effective parameters. Using the pore-scale distribution parameters as input, diffusion control equations are established for nanopores and micropores on the pore-scale cells, and periodic boundary conditions are applied. The Knudsen diffusion coefficient for nanopores and the Fick diffusion coefficient for micropores are obtained through numerical solution. After being transferred to the Darcy-scale equations through first-order asymptotic expansion, they together with the fracture convection mass transfer coefficients to form a multi-scale effective transport coefficient matrix.

[0013] The Knudsen diffusion coefficient of the nanopores was obtained by preparing a core sample with a diameter of 25 mm and a length of 50 mm, conducting a steady-state seepage experiment at 25°C using nitrogen as the medium, and covering the pore pressure gradient. ~ At least 8 sets of flow rate data points were collected within the MPa / m range. The linear relationship between the inverse of pressure and apparent permeability was fitted using the Klinkenberg method. The intrinsic permeability was obtained by extrapolation to infinite pressure. The Knudsen diffusion coefficient of the nanopores was then calculated by combining the difference between the intrinsic permeability and the apparent permeability with the pore size distribution parameters.

[0014] The discrete fracture network random generation technique refers to using the Monte Carlo method to independently sample according to the multi-parameter probability distribution function of fracture trace length, fracture aperture, fracture dip direction, and fracture inclination angle, and randomly generate a geometric set of fractures that satisfy statistical laws within the three-dimensional reservoir volume. Isolated fracture units are then eliminated through graph theory connectivity analysis, retaining the fracture network that constitutes the seepage channel. The multi-parameter probability distribution of fracture trace length, fracture aperture, fracture dip direction, and fracture inclination angle is jointly constrained by core fracture statistics, well logging imaging interpretation, and microseismic location data. After at least 200 independent Monte Carlo samplings, the convergence of sample statistical moments is used as the basis for determining the parameter distribution.

[0015] The integrated Kalman filter data assimilation method refers to approximating the error covariance in the Bayesian update with the integrated covariance matrix of the fracture parameter sample set. The residual between the dynamic pressure observation data and the forward model prediction value is mapped to the fracture aperture, fracture direction and fracture connectivity parameter space through the observation operator. Linear updates are performed synchronously for each integrated member. The iteration continues until the statistical distribution of the fracture parameter sample set converges, and the posterior probability distribution of fracture parameters, the equivalent permeability tensor and the matrix fracture quality exchange function are output.

[0016] Specifically, the adaptive multigrid enzyme kinetics operator splitting acceleration algorithm decomposes the matrix crack coupled diffusion reaction equation into diffusion sub-steps and enzyme-catalyzed reaction sub-steps in each time step according to the Strang symmetric splitting strategy. The diffusion sub-step constructs a 5-layer grid hierarchy using an algebraic multigrid framework. In the finest layer, Gauss-Seidel iteration is used to eliminate high-frequency error components, and in the coarsest layer, a direct solver is used to correct low-frequency error components. The residuals and corrections are transferred between layers using standard constraint operators and extension operators. The enzyme-catalyzed reaction sub-step uses a Krylov subspace exponential time integrator to process the rigid nonlinear source terms of the Michaelis-Menten equations for each component of the ESS-100 and ELS bioenzymes.

[0017] The dimensionless form of the Michaelis-Menten equation is as follows: ;in For the rate of enzyme-catalyzed reaction, For the maximum reaction rate, Substrate concentration, The Michaelis constant is used to determine the maximum reaction rates and Michaelis constants of each component of the ESS-100 bioenzyme and the ELS bioenzyme. These values ​​were obtained from isothermal and isobaric core enzymatic reaction experiments at reservoir temperatures of 60–90 °C and formation water salinity of 10–200 g / L. The results were obtained by nonlinear least squares fitting of the substrate concentration versus the enzyme reaction rate curve.

[0018] The enzyme concentration spatiotemporal evolution prediction model is a multi-scale spatiotemporal decomposition convolution-Transformer hybrid network. The model input is a time-series snapshot sequence on a three-dimensional reservoir grid, with the multi-scale effective transmission coefficient matrix, equivalent permeability tensor, and numerical solutions of the grid node concentration field as channel features. The spatial dimension is designed with three levels of parallel convolution branches. The output feature maps of the three parallel convolution branches are uniformly upsampled to an intermediate resolution by a three-dimensional cross-scale feature pyramid fusion module and then spliced ​​along the channel dimension. The spliced ​​features are input along the time axis to a causal Transformer encoder. The encoder output is restored to the three-dimensional concentration field temporal prediction result by a fully connected decoding layer.

[0019] The steps for establishing the training dataset for the enzyme concentration spatiotemporal evolution prediction model specifically include: using an adaptive multigrid enzyme kinetic operator splitting acceleration algorithm to train the dataset for a coverage porosity of 5%–25% and a crack density of 0.1–2.0. 300 parameter combinations of ESS-100 and ELS bioenzymes with concentrations ranging from 10 to 5000 mg / L were run until convergence. A complete three-dimensional concentration field time-series snapshot of each combination was saved as a labeled sample. The training, validation, and test sets were divided in an 8:1:1 ratio. The training steps specifically included: using the Adam optimizer at an initial learning rate... The training process is repeated for 300 epochs. If the validation set loss does not decrease after every 50 epochs, the learning rate is reduced to 0.5 times the current value.

[0020] The formula for calculating the value of the concentration field adaptive control function is as follows: ;in This represents the value of the concentration field adaptive control function. To predict residuals, The posterior standard deviation of the crack parameters. For coupling residuals, , , These are the corresponding reference standard values. , , These are the weighting coefficients.

[0021] Specifically, the physical information loss weight coefficient is dynamically adjusted according to the interval to which the concentration field adaptive control function value belongs. When the physical information loss weight coefficient is increased to 2.0 times its current value, a self-correcting iterative jump mechanism is triggered to recalculate the current time step; when At that time, the physical information loss weighting coefficient remains unchanged, and the number of local mesh layers at the crack matrix interface is increased by 1 layer only; when At that time, the physical information loss weight coefficient is reduced to 0.5 times the current value and the mesh is coarsened for uniform regions.

[0022] The reference standard values ​​for the predicted residuals, the posterior standard deviations of fracture parameters, the reference standard values ​​for the coupled residuals, and the weighting coefficients were obtained as follows: using 30 sets of core displacement experimental data covering different porosities and fracture densities as a benchmark, single-factor sensitivity scanning experiments were conducted on the predicted residuals, the posterior standard deviations of fracture parameters, and the coupled residuals, respectively. The contribution ratio of each factor to the mean square error of the final three-dimensional concentration field time-series prediction results was used to determine the weighting coefficients. , , The statistical mean of each input variable in the convergent state is used as the corresponding reference standard value.

[0023] Among them, the ESS-100 bioenzyme is a complex system composed of protein, alcohol dehydrogenase and xylanase, which migrates in the dense reservoir nanopores mainly through nanopore Knudsen diffusion; the ELS bioenzyme is a complex system composed of protease and lipase, in which the protease and lipase components migrate with the fracture convection mass transfer in the fracture channel and undergo enzymatic reactions on the fracture wall; the alcohol dehydrogenase and xylanase components in the ESS-100 bioenzyme, and the protease and lipase components in the ELS bioenzyme, are parameterized by their respective maximum reaction rates and Michaelis constants in the Michaelis-Menten equations of the enzymatic reaction substeps.

[0024] The equivalent permeability tensor describes the average seepage capacity along different spatial directions in a matrix fractured dual-pore medium in tensor form. Each component is obtained by homogenization calculation based on fracture aperture, fracture orientation, and fracture connectivity. The off-diagonal component of the tensor reflects the coupling effect of fracture anisotropy on seepage. The matrix fracture mass exchange function is a source-sink term function describing the molecular mass transfer of ESS-100 bio-enzymes and ELS bio-enzymes per unit volume of matrix and fracture per unit time. Its form is derived from the Warren-Root dual-pore model and parameterized by the fracture shape factor and the matrix effective diffusion coefficient component in the multi-scale effective transport coefficient matrix.

[0025] The causal mask is an upper triangular masking matrix that resets the position weights after the current moment to zero in the multi-head self-attention calculation of the causal Transformer encoder; the Strang symmetric splitting is the alternation of the diffusion sub-step and the enzymatic reaction sub-step in the operator splitting in a symmetric order of half-step-full-step-half-step; the posterior standard deviation of the crack parameters is the standard deviation of the sample set of crack aperture, crack direction and crack connectivity parameters after the iterative convergence of the integrated Kalman filter data assimilation method.

[0026] This invention combines multi-scale homogenization theory, discrete fracture network random generation technology, and integrated Kalman filter data assimilation method to construct a multi-scale effective transport coefficient matrix from nanopores to fracture scale. It replaces analytical covariance propagation with integrated statistics, mathematically eliminating the root cause of ill-conditioned Jacobian matrix in high-dimensional fracture parameter inversion. Furthermore, this invention employs an adaptive multi-grid enzyme kinetic operator splitting acceleration algorithm to decouple diffusion and enzymatic reactions, removing the constraint of rigid coupling on the time step and enabling efficient convergence of the coupled equations while concentrating computational resources in regions of abrupt concentration gradient changes. Simultaneously, a multi-scale spatiotemporal decomposition convolution-Transformer hybrid network embeds the diffusion equation residuals into the training objective, ensuring that the three-dimensional concentration field time-series prediction results satisfy physical laws under sparse labeling conditions. In summary, this invention solves the technical problems mentioned in the background art, such as the ill-conditioned equations and insufficient physical self-consistency of cross-scale concentration field time-series predictions in the multi-scale mass transfer coupling calculation of biological enzymes in the dual-porosity medium of tight reservoir matrix fractures. Attached Figure Description

[0027] Figure 1 This is a flowchart of the method of the present invention.

[0028] Figure 2 The graph shows the variation of the posterior standard deviation of crack opening with the number of iterations during the integrated Kalman filter data assimilation iteration process.

[0029] Figure 3 Three-dimensional concentration field distribution of ESS-100 bio-enzyme ethanol dehydrogenase component and ELS bio-enzyme protease component at different injection times in a matrix fracture dual-pore medium. Detailed Implementation

[0030] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below.

[0031] like Figure 1 The diagram shown is a flowchart of a method for establishing a biological enzyme diffusion model of a tight reservoir matrix fracture, provided by the present invention. This method includes the following steps:

[0032] S01. Collect core samples from tight reservoirs and obtain pore-scale distribution parameters from nanopores to millimeter-scale fractures through nitrogen adsorption experiments and mercury intrusion porosimetry experiments. Based on the multi-scale homogenization theory, the Knudsen diffusion coefficient of nanopores, the Fick diffusion coefficient of micropores, and the fracture convection mass transfer coefficient are asymptotically expanded and uploaded to the Darcy scale to construct a multi-scale effective transport coefficient matrix of hierarchical transfer.

[0033] S02. Using the discrete fracture network random generation technique, the pore scale distribution parameters obtained in S01 are used to construct the fracture geometric topology integrated sample. The dynamic pressure observation data is injected into the Monte Carlo forward propagation process through the integrated Kalman filter data assimilation method to constrain the posterior probability distribution of fracture aperture, fracture orientation and fracture connectivity, and output the equivalent permeability tensor and matrix fracture mass exchange function.

[0034] S03. Using the equivalent permeability tensor and matrix fracture mass exchange function obtained in S02 as input, the matrix fracture coupled diffusion reaction equation is decomposed into diffusion sub-steps and enzyme-catalyzed reaction sub-steps. The adaptive multigrid enzyme kinetic operator splitting acceleration algorithm is used to solve the problem iteratively on a multi-grid. After each complete operator splitting cycle, the coupling residual is calculated and the number of local grid layers and the number of smoothing iterations are dynamically adjusted. The numerical solution of the grid node concentration field at each time step is output.

[0035] S04. Input the multi-scale effective transmission coefficient matrix obtained in S01, the equivalent permeability tensor and matrix fracture mass exchange function obtained in S02, and the numerical solution of the grid node concentration field obtained in S03 into the enzyme concentration spatiotemporal evolution prediction model. The enzyme concentration spatiotemporal evolution prediction model outputs the three-dimensional concentration field time-series prediction results of ESS-100 bioenzyme and ELS bioenzyme in the matrix fracture dual-porosity medium.

[0036] S05. Using the coupling residual obtained in S03, the posterior standard deviation of the crack parameters obtained in S02, and the prediction residual of the three-dimensional concentration field time-series prediction result obtained in S04 as inputs, calculate the concentration field adaptive control function value. According to the interval to which the concentration field adaptive control function value belongs, dynamically adjust the physical information loss weight coefficient of the enzyme concentration spatiotemporal evolution prediction model, and re-drive the adaptive multigrid enzyme kinetic operator splitting acceleration algorithm of S03 to iterate until the global error tolerance is met.

[0037] S06. Based on the convergent three-dimensional concentration field time-series prediction results that meet the global error tolerance output in S05, and combined with the statistical results of the posterior probability distribution of the crack parameters obtained in S02, the spatial distribution of each component of ESS-100 bio-enzyme and ELS bio-enzyme in the matrix and cracks is evaluated, and the final matrix crack bio-enzyme diffusion model parameter set is output.

[0038] The multi-scale homogenization theory refers to a theoretical framework that uses asymptotic expansion as a mathematical tool to eliminate the microscopic heterogeneity in the pore-scale control equations through periodic cell integration, thereby obtaining macroscopic effective parameters. In this scheme, its specific implementation is as follows: using the pore-scale distribution parameters obtained from core mercury intrusion porosimetry as input, diffusion control equations are established for nanopores and micropores respectively on the pore-scale cells, and periodic boundary conditions are applied. The Knudsen diffusion coefficient for nanopores and the Fick diffusion coefficient for micropores are obtained through numerical solution. Then, through first-order asymptotic expansion, these coefficients are transferred to the Darcy-scale equations, forming a multi-scale effective transport coefficient matrix together with the fracture convection mass transfer coefficient. The technical effect of the multi-scale homogenization theory is that it eliminates the mass transfer error when a single-scale model cannot capture the coexistence of the nanopore Knudsen effect and the fracture convection effect through rigorous mathematical scale separation. This ensures that the control equations for different pore scales are coordinated within a unified mathematical framework, guaranteeing the physical consistency of the matrix fracture bioenzyme diffusion model in multi-scale media.

[0039] The Knudsen diffusion coefficient of the nanopores was obtained by preparing a core sample with a diameter of 25 mm and a length of 50 mm, conducting a steady-state seepage experiment at 25°C using nitrogen as the medium, and covering the pore pressure gradient. ~ At least 8 sets of flow data points were collected within the MPa / m range. The linear relationship between the inverse of pressure and apparent permeability was fitted using the Klinkenberg method. The intrinsic permeability was obtained by extrapolation to infinite pressure. The Knudsen diffusion coefficient of the nanopores was then calculated by combining the difference between the intrinsic and apparent permeability with the pore size distribution parameters. The asymptotic expansion refers to expressing the physical quantity as a series with the ratio of pore size to macroscopic size as a small parameter. At each expansion order, the corresponding periodic cell equation was established and solved. The results of the first two orders were used to form the components in the multi-scale effective transport coefficient matrix.

[0040] The discrete fracture network random generation technique refers to using the Monte Carlo method to independently sample multiple parameters (fracture trace length, fracture aperture, fracture dip direction, and fracture inclination angle) according to their probability distribution functions. This generates a set of geometrically consistent fractures within the three-dimensional reservoir volume, and isolated fracture units are eliminated through graph connectivity analysis, retaining the fracture network that forms the seepage channels. The probability distribution of these parameters is jointly constrained by core fracture statistics, well logging interpretation, and microseismic location data. After at least 200 independent Monte Carlo samplings, the convergence of the sample statistical moments serves as the basis for determining the parameter distribution. The technical benefits of this discrete fracture network random generation technique are: by replacing the deterministic fracture geometry description with a probability distribution, the uncertainty of the fracture network is explicitly expressed in the form of a sample set, providing prior sample support for the integrated Kalman filter data assimilation method, and fundamentally avoiding the ill-conditioned Jacobian matrix problem caused by a single deterministic fracture model.

[0041] The integrated Kalman filter data assimilation method refers to approximating the error covariance in Bayesian updates with the integrated covariance matrix of the fracture parameter sample set. Through observation operators, the residuals between dynamic pressure observation data and forward model predictions are mapped to the parameter space of fracture aperture, fracture orientation, and fracture connectivity. Linear updates are synchronously performed on each integrated member, iterating until the statistical distribution of the fracture parameter sample set converges. The resulting output includes the posterior probability distribution of fracture parameters, the equivalent permeability tensor, and the matrix-fracture mass exchange function. The technical benefits of this integrated Kalman filter data assimilation method are: replacing analytical covariance propagation with integrated statistics, maintaining computational feasibility in high-dimensional ill-posed problems where the dimension of fracture parameters is much larger than the amount of dynamic pressure observation data; and continuously compressing the posterior uncertainty of fracture parameters using dynamic pressure observation data, gradually approximating the estimated value of the equivalent permeability tensor to the true state of the reservoir, thus completely eliminating the ill-conditioned nature of the equations caused by fracture topological uncertainties.

[0042] The specific implementation of the adaptive multigrid enzyme kinetics operator splitting acceleration algorithm is as follows: The matrix crack coupled diffusion reaction equation is decomposed into diffusion sub-steps and enzyme-catalyzed reaction sub-steps alternately in each time step according to the Strang symmetric splitting strategy. The diffusion sub-step constructs a 5-layer grid hierarchy from the finest to the coarsest grid using an algebraic multigrid framework. In the finest layer, Gauss-Seidel iteration is used to eliminate high-frequency error components, and in the coarsest layer, a direct solver is used to correct low-frequency error components. Residuals and corrections are propagated between layers using standard constraint operators and extension operators. The enzyme-catalyzed reaction sub-step uses a Krylov subspace exponential time integrator to process the rigid nonlinear source terms of the Michaelis-Menten equations for each component of the ESS-100 and ELS bioenzymes. After each complete operator splitting cycle, the algorithm is further refined using... The norm is used to calculate the spatial distribution of the coupling residuals. At the fracture-matrix interface, the number of grid layers is locally increased for grid cells where the coupling residuals exceed a preset coupling residual threshold. In uniform regions where the coupling residuals are below the preset coupling residual threshold, the number of smoothing iterations is reduced. The preset coupling residual threshold is obtained through iterative experiments using 10 sets of core simulation examples with different porosities and fracture densities. The 95th percentile of the relative error of the concentration field between coarse and fine grids for each example is used as the basis for setting the preset coupling residual threshold. The adaptive multigrid enzyme kinetic operator splitting acceleration algorithm brings the following technical effects to the scheme: Strang symmetric splitting allows diffusion and enzymatic reactions to be processed independently under their respective optimal numerical formats, eliminating the time step limitation imposed by rigid coupling; the adaptive grid hierarchy strategy concentrates computational resources in the fracture-matrix interface region where the concentration gradient changes sharply, and automatically coarsens in uniform regions, minimizing the overall floating-point computation while meeting the global error tolerance.

[0043] The dimensionless form of the Michaelis-Menten equation is as follows: ;in For the rate of enzyme-catalyzed reaction ( ), For the maximum reaction rate ( ), substrate concentration ( ), Michaelis constant ( The maximum reaction rates and Michaelis constants of the ESS-100 bioenzyme and each component of the ELS bioenzyme were determined by isothermal and isobaric core enzymatic reaction experiments in the reservoir temperature range of 60-90℃ and formation water salinity range of 10-200 g / L. The results were obtained by nonlinear least square fitting of substrate concentration and enzyme reaction rate curves.

[0044] The enzyme concentration spatiotemporal evolution prediction model is a multi-scale spatiotemporal decomposition convolution-Transformer hybrid network. The specific structure of the model is as follows: the model input is a time-series snapshot sequence of channel features on a three-dimensional reservoir grid, with the multi-scale effective transport coefficient matrix, equivalent permeability tensor, and numerical solutions of grid node concentration fields as channel features. A three-level parallel convolution branch is designed in the spatial dimension. The first branch uses deformable convolution kernels with receptive fields corresponding to the nanopore scale to extract local pore geometric features of nanopores. The deformable convolution adds a learnable displacement field to the sampling offset of the standard convolution, adaptively sensing irregular pore deformation. The second branch uses deformable convolution kernels with receptive fields corresponding to the micropore scale to extract micropore network features. The third branch uses deformable convolution kernels with receptive fields corresponding to the fracture scale to extract fracture geometric features. The output of the third branch is expanded in the channel dimension to accommodate fracture anisotropy information. The output feature maps of the three parallel convolution branches are uniformly upsampled to an intermediate resolution by a three-dimensional cross-scale feature pyramid fusion module and then stitched along the channel dimension. The three-dimensional cross-scale feature pyramid fusion module completes spatial alignment using three-dimensional bilinear interpolation. Convolution completes channel compression; the concatenated features are input along the time axis into a causal Transformer encoder. This encoder uses a causal mask to ensure that the prediction time only focuses on historical information. The position encoding is constructed by concatenating the normalized value of the injection time with the historical sequence of formation pressure. A multi-head self-attention mechanism captures long-range dependencies across time steps. The encoder output is then restored to the three-dimensional concentration field time-series prediction result through a fully connected decoding layer. The model incorporates a self-correcting iterative jump mechanism, calculating the prediction residual after each time step output. When the prediction residual exceeds the self-correction threshold, the model backtracks to... The self-calibration threshold was determined by cross-validation of 50 sets of labeled laboratory core flow experimental data at the previous time step. The total loss function was constructed by summing the diffusion equation residuals, which were weighted by the mean square error of the data fitting and the weighting coefficient of the physical information loss. The diffusion equation residuals incorporated the residuals of the diffusion equation at the collocation points into backpropagation, constraining the model to satisfy physical laws under sparse labeling conditions. The steps for establishing the training dataset for the enzyme concentration spatiotemporal evolution prediction model specifically included: using an adaptive multigrid enzyme kinetic operator splitting acceleration algorithm to train the model with a porosity of 5%–25% and a fracture density of 0.1–2.0. The ESS-100 and ELS bioenzymes, with concentrations ranging from 10 to 5000 mg / L, were run to convergence with 300 parameter combinations. A complete three-dimensional concentration field time-series snapshot of each combination was saved as a labeled sample. The training, validation, and test sets were divided in an 8:1:1 ratio. The training steps for the enzyme concentration spatiotemporal evolution prediction model specifically included: using the Adam optimizer at an initial learning rate... The training process is repeated for 300 rounds. Every 50 rounds, if the validation set loss does not decrease, the learning rate is reduced to 0.5 times the current value. The final model parameters are the weights corresponding to the minimum sum of the diffusion equation residuals and the mean square error of the data fitting on the validation set. The technical effects of the enzyme concentration spatiotemporal evolution prediction model are as follows: the three-level parallel deformable convolutional branches synchronously encode the spatial heterogeneity of the three scales of nanopores, micropores and cracks, avoiding the underexpression of cross-scale pore structures by single-scale receptive fields; the causal Transformer encoder captures long-range temporal dependencies in the injection history with a multi-head self-attention mechanism, making up for the gradient vanishing problem of traditional recurrent networks on long sequences; the diffusion equation residuals weighted by the physical information loss weight coefficients embed the diffusion equation constraints into the training objective, so that the model still maintains physically consistent three-dimensional concentration field temporal prediction results in the sparse reservoir region of the training samples. Overall, the spatiotemporal prediction of the matrix crack biological enzyme diffusion model maintains physical self-consistency under multi-scale geometry and sparse observation conditions.

[0045] The formula for calculating the value of the concentration field adaptive control function is as follows: ;in The value of the concentration field adaptive control function (dimensionless) To predict residuals ( ), The posterior standard deviation (m) of the crack parameters. For coupling residuals ( ), To predict the reference standard value of the residual ( ), The reference standard value (m) is the posterior standard deviation of the crack parameters. The reference standard value for coupling residuals ( ), , , For weighting coefficients (dimensionless); when When the physical information loss weight coefficient is increased to 2.0 times its current value, a self-correcting iterative jump mechanism is triggered to recalculate the current time step; when At that time, the physical information loss weighting coefficient remains unchanged, and the number of local mesh layers at the crack matrix interface is increased by 1 layer only; when At that time, the physical information loss weighting coefficient is reduced to 0.5 times the current value, and the mesh is coarsened for uniform regions. The method for obtaining the prediction residual reference standard value, fracture parameter posterior standard deviation reference standard value, coupling residual reference standard value, and weighting coefficient is as follows: based on 30 sets of core displacement experimental data covering different porosities and fracture densities, single-factor sensitivity scanning experiments are performed on the prediction residual, fracture parameter posterior standard deviation, and coupling residual, respectively, and the contribution ratio of each factor to the mean square error of the final three-dimensional concentration field time series prediction result is used to determine the weighting coefficient. , , The statistical mean of each input variable in the convergent state is used as the corresponding reference standard value.

[0046] Among them, the ESS-100 bioenzyme is a composite system composed of protein, alcohol dehydrogenase, and xylanase. It migrates within the nanopores of dense reservoirs primarily via Knudsen diffusion. The alcohol dehydrogenase component acts on the hydroxyl functional groups of residual organic matter in the reservoir, while the xylanase component hydrolyzes xylan polysaccharides adsorbed on the pore walls. The alcohol dehydrogenase and xylanase components synergistically reduce interfacial tension in the matrix pores, thereby improving the effective permeability of the pores. The ELS bioenzyme is a composite system composed of protease and lipase. The protease component hydrolyzes and blocks... Protein-like deposits in the pore throat; lipase components catalyze the cleavage of long-chain fatty acid ester bonds adsorbed on the reservoir rock surface; protease and lipase components migrate with fracture convection mass transfer in the fracture channel and undergo enzymatic reactions on the fracture wall, improving the wetting state by reducing the hydrophobicity of the fracture surface; the alcohol dehydrogenase and xylanase components in ESS-100 bioenzyme, and the protease and lipase components in ELS bioenzyme, are parameterized by their respective maximum reaction rates and Michaelis constants in the Michaelis-Menten equations of the enzymatic reaction substeps.

[0047] The equivalent permeability tensor refers to the average seepage capacity along different spatial directions in a matrix fractured dual-pore medium, described in tensor form. Its components are obtained through homogenization calculations based on fracture aperture, fracture orientation, and fracture connectivity. The off-diagonal components of the tensor reflect the coupling effect of fracture anisotropy on seepage. The matrix fracture mass exchange function describes the source-sink term function describing the molecular mass transfer between ESS-100 and ELS bioenzymes per unit volume of matrix and fracture per unit time. Its form is derived from the Warren-Root dual-pore model and parameterized by the fracture shape factor and the matrix effective diffusion coefficient component in the multi-scale effective transport coefficient matrix. The three-dimensional cross-scale feature pyramid fusion module refers to aligning three-dimensional feature maps of different resolutions through three-dimensional bilinear interpolation space and then stitching them together in the channel dimension. The multi-scale feature fusion structure of convolutional compression transmits high-level semantic features to low-level spatial detail features via a top-down path; the causal mask refers to the upper triangular masking matrix that resets the position weights after the current time step to zero in the multi-head self-attention calculation of the causal Transformer encoder, ensuring the temporal causality of information flow during prediction; the Strang symmetric splitting refers to alternating the diffusion sub-step and the enzymatic reaction sub-step in the operator splitting in a symmetric order of half-step-full-step-half-step, thereby increasing the order of time integral error to second order; the posterior standard deviation of the crack parameters refers to the standard deviation of the sample set of crack aperture, crack direction, and crack connectivity parameters after the iterative convergence of the integrated Kalman filter data assimilation method, reflecting the degree of dispersion of the posterior probability distribution of crack parameters.

[0048] The specific implementation of step S01 is as follows: Technicians retrieve representative core samples with a diameter of 25 mm and a length of 50 mm from the target tight reservoir section. First, a low-temperature liquid nitrogen adsorption experiment is used to obtain the pore size distribution of the nanopore to micropore range (2–200 nm). Then, a mercury intrusion porosimetry experiment is used to obtain the pore size distribution of the micropore to millimeter-scale fracture range (100 nm–1 mm). The two sets of distribution data are interpolated and stitched together within the overlapping pore size range to form a set of pore-scale distribution parameters covering the entire scale. In the nanopore Knudsen diffusion coefficient acquisition stage, nitrogen gas is used as the seepage medium, and different pore pressure gradients are applied under a constant temperature of 25°C. ~ At least eight sets of steady-state flow rate data points were collected (MPa / m). The apparent permeability was plotted against the reciprocal of pressure using the Klinkenberg linear fitting method. The intrinsic permeability was obtained by extrapolating to where the reciprocal of pressure approaches zero. The Knudsen diffusion coefficient for nanopores was then calculated using the Knudsen number in conjunction with pore-scale distribution parameters. Subsequently, steady-state diffusion control equations were established for nanopores and micropores at the pore-scale cell level, respectively. Three-dimensional periodic boundary conditions were applied, and the cell problems were numerically solved using the finite element method. The results were then integrated using first-order asymptotic expansion to obtain the effective diffusion coefficients for the corresponding scales. These effective diffusion coefficients were then assembled with the fracture convection mass transfer coefficients identified from pressure transient tests to form a multi-scale effective transport coefficient matrix. The purpose of this step is to provide a physically consistent parameter basis for subsequent mass transfer calculations at all scales, ensuring that the nanopore Knudsen effect and the fracture convection effect are correctly described within a unified mathematical framework.

[0049] The specific implementation of step S02 is as follows: Using the pore-scale distribution parameters obtained in step S01 as prior constraints, marginal probability distribution functions for fracture trace length, fracture aperture, fracture dip direction, and fracture inclination are extracted from core fracture statistics, well logging imaging interpretation, and microseismic location data. The Monte Carlo method is used to independently sample these four parameters, generating a set of fracture geometry that satisfies statistical laws within the three-dimensional reservoir volume. The connectivity of fracture units is analyzed using graph theory adjacency matrices, eliminating isolated fractures not connected to the main seepage network and retaining effective seepage channels. Each independent sampling generates one ensemble member, repeated at least 200 times until the sample statistical moments (mean and variance) converge, forming the prior sample set for ensemble Kalman filtering. In the data assimilation stage, a forward flow simulator is used to calculate the predicted pressure of each ensemble member at the dynamic pressure observation time. The observation operator maps the residual between the measured pressure and the predicted pressure to the fracture aperture, fracture direction, and fracture connectivity parameter space through the ensemble covariance matrix. A linear Kalman update is performed on each ensemble member, iterating cyclically until the statistical distribution of fracture parameters for all ensemble members converges. After convergence, the equivalent permeability tensor components are calculated using the integrated mean. The matrix fracture mass exchange function is constructed using the Warren-Root dual-pore model combined with the fracture shape factor and the effective diffusion coefficient of the matrix. At the same time, the posterior standard deviation of the fracture parameters of each integrated member is saved for use in step S05.

[0050] The specific implementation of step S03 is as follows: Substitute the equivalent permeability tensor and matrix fracture mass exchange function output from step S02 into the matrix fracture coupled diffusion reaction equation. Using the Strang symmetric splitting strategy, each time step is decomposed into a symmetric sequence of half-step diffusion sub-steps, full-step enzymatic reaction sub-steps, and half-step diffusion sub-steps. The diffusion sub-step constructs a five-layer hierarchy from the finest grid (grid spacing corresponding to fracture aperture magnitude) to the coarsest grid within an algebraic multigrid framework. The finest layer uses Gauss-Seidel iteration to smooth high-frequency errors, and layer-by-layer, the residuals are propagated to the coarsest layer using constraint operators. In the coarsest layer, low-frequency errors are corrected using a direct solver, and then the correction is propagated back to the finest layer layer by layer using an extension operator, completing one V-cycle. The enzymatic reaction sub-step, targeting the alcohol dehydrogenase component and xylanase component in ESS-100 bioenzymes, and the protease component and lipase component in ELS bioenzymes, is substituted into the corresponding Michaelis-Menten equations. A Krylov subspace exponential time integrator is used to process the rigid nonlinear source terms, with a reference step size of [missing information]. On the order of s. After each complete operator splitting loop, the value is... The norm was used to calculate the spatial distribution of the global coupling residuals. At the fracture-matrix interface, the coupling residuals exceeded a preset coupling residual threshold (determined by the 95th quantile of 10 sets of examples with different porosities and fracture densities; the reference value is approximately...). The number of grid layers is locally increased in grid cells (on the order of magnitude), the number of smoothing iterations is reduced in uniform regions where the coupling residual is below the threshold, computational resources are adaptively allocated, and the numerical solution of the grid node concentration field at each time step is iteratively output.

[0051] The specific implementation of step S04 is as follows: The outputs of steps S01 to S03, namely the multi-scale effective transmission coefficient matrix, equivalent permeability tensor, matrix fracture mass exchange function, and numerical solutions of grid node concentration fields for all time steps, are spliced ​​into a multi-channel time-series snapshot sequence on a three-dimensional reservoir grid, and input into the enzyme concentration spatiotemporal evolution prediction model. The three-level parallel convolutional branches of this model extract spatial pore geometric features at each scale using deformable convolutional kernels corresponding to the receptive field at the nanopore scale, micropore scale, and fracture scale, respectively. The outputs of the three branches are aligned to the intermediate resolution by a three-dimensional cross-scale feature pyramid fusion module using three-dimensional bilinear interpolation and then spliced ​​along the channel dimension. The fused time-series features are input along the time axis into a causal Transformer encoder. The position encoding is composed of the normalized value of the injection time and the formation pressure history sequence. A multi-head self-attention mechanism captures long-range dependencies across time steps. The encoder output is restored to the three-dimensional concentration field time-series prediction result by a fully connected decoding layer. The model's internal self-calibrating iterative jump mechanism calculates the prediction residual after each time step output. When the prediction residual exceeds the self-calibration threshold (determined by cross-validation of 50 sets of laboratory core flow experimental data, with a reference value of approximately...), the self-calibration threshold is reached. At times of magnitude (on the order of magnitude), backtracking and recalculating ensures the physical continuity of the output at each time step. The purpose of this step is to use artificial intelligence methods to efficiently replace pure numerical solutions for time series prediction under a large number of parameter combinations, while maintaining the physical consistency of the prediction results through physical information loss constraints.

[0052] The specific implementation of step S05 is as follows: the coupling residual output in step S03 is... The posterior standard deviation of the crack parameters output in step S02 The predicted residual output from step S04 Substitute into the concentration field adaptive control function Calculate the dimensionless control function value The reference standard values ​​and weighting coefficients were determined by single-factor sensitivity scanning of 30 sets of core displacement experimental data, with the reference value range being [range missing]. Approximately 0.4 Approximately 0.3 Approximately 0.3. According to... Differential regulation will be implemented within the relevant interval: when Increase the weighting coefficient for physical information loss to 2.0 times the current value and trigger self-correcting backtracking; when While keeping the weights constant, only the number of mesh layers at the interface is increased by 1; when The physical information loss weight coefficient is reduced to 0.5 times its current value, and the mesh in the uniform region is coarsened. After adjustment, the adaptive multigrid enzyme kinetic operator splitting acceleration algorithm in step S03 is re-driven iteratively until the global coupling residuals are all below the global error tolerance (reference value). (on a scale of magnitude), achieving closed-loop collaborative convergence between the solver and the prediction model.

[0053] The specific implementation of step S06 is as follows: Based on the three-dimensional concentration field time-series prediction results converged in step S05, and combined with the statistical results of the posterior probability distribution of fracture parameters saved in step S02, the spatial distribution of alcohol dehydrogenase and xylanase components in the matrix nanopores of ESS-100 bioenzymes, and the spatial distribution of protease and lipase components in the fracture channels of ELS bioenzymes are statistically evaluated. The mean and standard deviation fields of the concentration fields of each component are calculated, the concentration gradient distribution characteristics at the fracture matrix interface are extracted, and the migration front position and distribution uniformity of bioenzymes in media with different pore scales are determined. The components of the converged multi-scale effective transport coefficient matrix, the components of the equivalent permeability tensor, the matrix fracture mass exchange function parameters, the maximum reaction rate and Michaelis constant of each component in the Michaelis-Menten equation, and the convergence statistics of the three-dimensional concentration field time-series prediction results are uniformly encapsulated and the final matrix fracture bioenzyme diffusion model parameter set is output, providing a model basis for the design of subsequent reservoir production enhancement schemes.

[0054] It should be noted that the key technologies of this invention include: the multi-scale homogenization theory eliminates pore-scale heterogeneity within periodic cell integrals through asymptotic expansion, enabling the Knudsen effect of nanopores and the convection effect of fractures to be rigorously expressed within the same mathematical framework, thereby eliminating the systematic errors generated by single-scale models in the cross-scale mass transfer coefficient transfer; the combination of discrete fracture network random generation technology and integrated Kalman filter data assimilation method, using the sample set covariance matrix instead of analytical covariance propagation, continuously compresses the posterior distribution of fracture parameters from dynamic pressure observation data, completely avoiding the Jacobian error caused by deterministic fracture models in high-dimensional parameter inversion. The algorithm addresses the problem of ill-conditioned matrices. An adaptive multigrid enzyme kinetics operator splitting acceleration algorithm decouples diffusion and enzymatic reactions through Strang-symmetric splitting, eliminates cross-frequency errors using an algebraic multigrid framework, and adaptively concentrates computational resources on the fracture matrix interface region with abruptly changing concentration gradients. An adaptive concentration field control function quantifies the error indices of the three processes, driving the solver and AI prediction model to converge iteratively. This forms a closed-loop collaborative mechanism among multi-scale mass transfer calculation, uncertainty quantification, and physical constraint prediction, achieving a physically self-consistent and efficient solution under conditions of coexisting fracture topological uncertainty and multi-scale pore structure.

[0055] It should be noted that the present invention also solves the following technical problems: In the modeling of bio-enzyme diffusion in tight reservoirs, ESS-100 bio-enzymes and ELS bio-enzymes have significantly different enzyme kinetic parameters, and the dominant mass transfer mechanisms of each component in the nanoporous matrix and fracture channels are essentially different. Traditional single-component diffusion models cannot simultaneously describe the coupled reaction migration behavior of multi-component enzymes in heterogeneous media, resulting in mutual interference of the time-series prediction results of the concentration fields of each component and a lack of component resolution. This invention parameterizes the alcohol dehydrogenase and xylanase components of ESS-100 bioenzymes, and the protease and lipase components of ELS bioenzymes, respectively, with their maximum reaction rates and Michaelis constants in the enzymatic reaction substeps of the Michaelis-Menten equation. Furthermore, it synchronously encodes the spatial heterogeneity of the nanoporous matrix and cracks in a three-level parallel convolution branch of a multi-scale spatiotemporal decomposition convolution-Transformer hybrid network. This achieves high-precision temporal prediction of coupled mass transfer and concentration fields of each component in multi-component bioenzymes within a matrix-crack dual-porous medium, solving the technical problem that the coupled reaction migration behavior of multi-component bioenzymes in heterogeneous dual-porous media cannot be uniformly modeled.

[0056] Specifically, the principle of this invention is as follows: The fundamental reason why this invention can solve the above-mentioned technical problems lies in the interconnected mathematical structure design at three levels. At the scale transfer level, the asymptotic expansion method eliminates the microscopic heterogeneity at the pore scale through periodic cell integration, enabling the Knudsen diffusion coefficient of nanopores and the Fick diffusion coefficient of micropores to be transferred to the Darcy scale through rigorous mathematical derivation rather than empirical fitting, thus ensuring the physical consistency of multi-scale mass transfer coefficients from the source. At the fracture uncertainty handling level, the integrated Kalman filter approximates the Bayesian error propagation with the sample set covariance matrix, continuously compressing the posterior distribution of fracture parameters from dynamic pressure observation data, so that the equivalent permeability tensor gradually converges to the true state of the reservoir, completely avoiding the ill-conditioned equations caused by deterministic fracture models. At the temporal prediction level, the causal Transformer encoder captures long-range temporal dependencies with a multi-head self-attention mechanism. The physical information loss weight coefficient embeds the diffusion equation constraint into backpropagation, enabling the model to produce physically consistent concentration field prediction results even in sparse regions of training samples. The concentration field adaptive control function further unifies and quantifies the coupling residual, crack parameter uncertainty and prediction residual, dynamically driving the operator split solver and the neural network prediction model to iterate collaboratively to the global error tolerance, forming a closed-loop adaptive correction mechanism.

[0057] The following provides a specific embodiment 1 of the present invention, and the specific implementation of each step in this embodiment 1 is described in detail below.

[0058] The specific implementation of step S01 is as follows: After collecting core samples from tight reservoirs, nitrogen adsorption experiments and mercury intrusion porosimetry experiments are carried out sequentially to obtain pore-scale distribution parameters from nanopores to millimeter-scale fractures. The steps for obtaining the Knudsen diffusion coefficient of nanopores are as follows: Step 1, prepare core samples with a diameter of 25 mm and a length of 50 mm; Step 2, conduct steady-state seepage experiments at 25℃ using nitrogen as the medium, under pore pressure gradient coverage. ~ Collect at least 8 sets of flow data points within the MPa / m range; Step 3: Fit the pressure reciprocal using the Klinkenberg method. With apparent penetration The linear relationship is expressed by the following formula:

[0059] ;

[0060] In the formula, Apparent permeability ( ) Intrinsic permeability ( ), by extrapolating to infinite pressure, that is get The Klinkenberg slip coefficient (Pa) is determined by the slope of the linear fit. The pore pressure is (Pa). Step 4, from... and The difference is combined with the pore size distribution parameters to back-calculate the Knudsen diffusion coefficient of the nanopores. ( Based on the multi-scale homogenization theory, diffusion control equations are established for nanopores and micropores at the pore-scale cell level, and periodic boundary conditions are applied. The Knudsen diffusion coefficient tensor of the nanopore is then expanded using a first-order asymptotic expansion. Fick diffusion coefficient tensor of micron pores Vector of mass transfer coefficients in the crack Upload to the Darcy scale to construct a multi-scale effective transport coefficient matrix. Nanopore Knudsen diffusion coefficient tensor. ( The formula for ) is expressed as follows:

[0061] ;

[0062] Micron-sized Fick diffusion coefficient tensor ( The formula for ) is expressed as follows:

[0063] ;

[0064] Crack convection mass transfer coefficient vector ( The formula for ) is expressed as follows:

[0065] ;

[0066] In the formula, For the tensor components of the Knudsen diffusion coefficient of nanopores ( (Originally obtained from the Klinkenberg experiment) For the tensor component of the Fick diffusion coefficient of a micrometer pore ( The pore size distribution parameters from the nitrogen adsorption experiment were obtained by combining the periodic cell values. The mass transfer coefficient of the crack convection in the first... Components of direction ( ), obtained by homogenizing the crack aperture and velocity distribution through integral; subscript Represents the spatial direction component. The multi-scale effective transmission coefficient matrix is ​​composed of... , and The set of graded transfer parameters, which together constitute the set, is called by each component according to the corresponding pore size in subsequent steps.

[0067] The specific implementation of step S02 is as follows: using the discrete crack network random generation technique, the crack trace length is determined by the Monte Carlo method. (m), crack aperture (m) Crack tendency (°), crack inclination angle The multi-parameter probability distribution function of (°) is independently sampled. After at least 200 independent samplings and using sample statistical moment convergence as the criterion, an integrated sample of the crack geometric topology is constructed. Subsequently, the dynamic pressure observation data is injected into the Monte Carlo forward propagation process using an integrated Kalman filter data assimilation method, and a linear update is performed on the crack parameter sample set. The update formula is expressed as follows:

[0068] ;

[0069] In the formula, For the first The first integrated member The crack parameter vector of the next iteration includes the crack aperture. , direction of cracks Connectivity with cracks After adding observation noise, the first The dynamic pressure observation vector (Pa) corresponding to each integrated member. The observation operator matrix maps the fracture parameter space to the pressure observation space, and is a dimensionless linear mapping matrix. For the first The Kalman gain matrix (dimensionless) for the next iteration is expressed as follows:

[0070] ;

[0071] In the formula, For the first The cross-covariance matrix of fracture parameters and predicted pressure at the next iteration has a row dimension representing the fracture parameter dimension and a column dimension representing the pressure observation dimension. For the first The ensemble covariance matrix for predicting pressure in the next iteration is a square matrix with dimensions equal to the dimension of the pressure observations. Let be the observation error covariance matrix, with dimensions and . The same characteristics reflect the statistical properties of measurement errors in dynamic pressure observation data. After iterating until the statistical distribution of the fracture parameter sample set converges, the equivalent permeability tensor is output. ( Mass exchange function with matrix cracks ( Equivalent permeability tensor The formula is expressed as follows:

[0072] ;

[0073] In the formula, For each component of the equivalent permeability tensor ( ), determined by crack aperture , direction of cracks Connectivity with cracks The off-diagonal components are obtained from the geometric parameters of the discrete crack network through homogenized integration after posterior estimation using the integrated Kalman filter data assimilation method. ( This reflects the coupling effect of fracture anisotropy on seepage. Matrix fracture mass exchange function. Derived from the Warren-Root two-hole model, the formula is expressed as follows:

[0074] ;

[0075] In the formula, Mass exchange function for matrix cracks ( ) Reference standard value for matrix crack mass exchange function ( ), used for dimensionless processing Crack shape factor ( Determined by crack geometry parameters The matrix effective diffusion coefficient component in the multi-scale effective transport coefficient matrix ( ), taking in the nanoporous region Corresponding components are taken in the micron-pore region. Corresponding components The concentration of biological enzymes in the crack ( ) The concentration of biological enzymes in the matrix ( ).

[0076] The specific implementation method of step S03 is: with and As input, the Strang symmetric splitting strategy was used to decompose the matrix fracture coupled diffusion reaction equation in a symmetric order of half-step diffusion, full-step enzymatic reaction, and half-step diffusion. The diffusion sub-step governing equation is expressed as follows:

[0077] ;

[0078] In the formula, Reference standard value for biological enzyme concentration ( ), used for dimensionless processing The time reference standard value (s) is used for dimensionless processing. The concentration of biological enzymes ( ) For time (s) gradient operator The effective diffusion tensor called according to the porosity region ( ), taking in the nanoporous region Take in the micron-pore region Mass exchange function for matrix cracks ( The result is obtained from step S02. The enzyme-catalyzed reaction substep uses the Michaelis-Menten equation to handle the rigid nonlinear source term, which is expressed in dimensionless form as follows:

[0079] ;

[0080] In the formula, For the rate of enzyme-catalyzed reaction ( ) For the maximum reaction rate ( ) substrate concentration ( ) Michaelis constant ( The data were obtained from isothermal and isobaric core enzymatic reaction experiments within a reservoir temperature range of 60–90℃ and formation water salinity of 10–200 g / L by nonlinear least squares fitting of substrate concentration versus enzymatic reaction rate curves. An adaptive multigrid framework was constructed with a 5-layer grid hierarchy. After each complete operator splitting cycle, the data was processed using... Norm calculation of coupled residuals Its formula is expressed as follows:

[0081] ;

[0082] In the formula, For coupling residuals ( ) For the computational domain ( ) The enzyme-catalyzed reaction rate calculated using the Michaelis-Menten equation ( ); the dimensions of each component of the integrand are all The square root of the integral over the volume of the computational domain has dimensions of . , Here, it is used as a measure of coupling imbalance in the sense of spatial average, with a reference standard value. The concentration field relative error between coarse and fine grids was statistically determined from 10 core simulation examples with different porosities and fracture densities. Dimensions are consistent. At the fracture matrix interface. For grid cells that exceed the preset coupling residual threshold, the number of grid layers is increased locally, and the number of smoothing iterations is reduced in uniform regions. The numerical solution of the grid node concentration field at each time step is output.

[0083] The specific implementation of step S04 is as follows: The multi-scale effective transmission coefficient matrix, , The numerical solution of the concentration field at the grid nodes is input into the spatiotemporal evolution prediction model of enzyme concentration. A three-level parallel deformable convolution branch extracts spatial features at the nanopore, micropore, and crack scales, respectively. After alignment by a three-dimensional cross-scale feature pyramid fusion module, these features are input along the time axis into a causal Transformer encoder to capture long-range dependencies across time steps. The final output is the three-dimensional concentration field time-series prediction results of ESS-100 and ELS bioenzymes in a matrix with dual porosity in cracks. The total loss function formula is expressed as follows:

[0084] ;

[0085] In the formula, The total loss function value ( ) The reference standard value for the loss function ( ), take the variance statistical mean of the numerical solutions of the concentration fields in the training set. For the first The model predicts the concentration at each point ( ) For the first The labeled concentration of each point ( ) The total number of points allocated to the training samples Number of physical points Fit dimensionless weighting coefficients to the data. The physical information loss weighting coefficient (dimensionless) For the first The enzyme reaction rate at each coordinate point, calculated using the Michaelis-Menten equation ( ),Depend on Calculation and acquisition, where and Take the value obtained from the core experiment fitting in step S03; the dimensions of the first data fitting term on the right side of the equation are... The second physical information term is the integral with dimensions of Its square dimension is To unify the two dimensions, Implicit in physical information items Normalization factor, in actual calculation The two values ​​are normalized separately and then weighted and summed.

[0086] The specific implementation of step S05 is as follows: using coupling residuals Posterior standard deviation of crack parameters Compared with the predicted residual Using the input, calculate the value of the adaptive control function for the concentration field. The formula is expressed as follows:

[0087] ;

[0088] In the formula, Values ​​of the concentration field adaptive control function (dimensionless) To predict residuals ( ), to predict concentration for the model Compared with the labeled concentration The spatial root mean square of the difference Let m be the posterior standard deviation of the crack parameters, and m be the crack aperture after iterative convergence of the integrated Kalman filter data assimilation method. , direction of cracks Connectivity with cracks The standard deviation of the parameter sample set reflects the degree of dispersion of the posterior probability distribution of crack parameters. For coupling residuals ( ) To predict the reference standard value of the residual ( ) Reference standard value (m) for the posterior standard deviation of crack parameters. The reference standard value for coupling residuals ( ) , , The weighting coefficients (dimensionless) were determined by single-factor sensitivity scanning of 30 sets of core displacement experimental data, and the default values ​​are satisfied. , , and , , The data was obtained by using 30 sets of core displacement experimental data covering different porosities and fracture densities as a benchmark. , and Single-factor sensitivity scanning experiments were conducted to determine the contribution ratio of each factor to the mean square error of the final three-dimensional concentration field time-series prediction results. , , The statistical mean of each input variable in the convergent state is used as the corresponding reference standard value. At that time, Increase to 2.0 times the current value and trigger a self-correcting iterative jump mechanism to recalculate the current time step; when At that time, keep The number of mesh layers remains unchanged, except that the number of mesh layers at the crack matrix interface is increased by one; when At that time, Reduce the value to 0.5 times the current value and coarse the mesh for the uniform region, then restart step S03 iteration until the global error tolerance is met.

[0089] The specific implementation of step S06 is as follows: Based on the convergent three-dimensional concentration field time-series prediction results that satisfy the global error tolerance output in step S05, and combined with the statistical results of the posterior probability distribution of crack parameters, the spatial distribution of alcohol dehydrogenase and xylanase components in ESS-100 bioenzymes, and protease and lipase components in ELS bioenzymes in the matrix and cracks is comprehensively evaluated, and the final matrix-crack bioenzyme diffusion model parameter set is output. This parameter set includes each component of the multi-scale effective transport coefficient matrix, , , and each component and The convergence estimate of .

[0090] To better understand and implement this invention, the following is a specific application scenario of the invention, Example 2: In order to verify the effect of the invention, the technicians set up a test environment and carried out a complete matrix fracture bio-enzyme diffusion modeling experiment by selecting a core sample of a tight sandstone reservoir. The applicability and computational convergence of the method of the invention under actual reservoir parameter conditions were systematically verified.

[0091] The core sample has a diameter of 25 mm and a length of 50 mm, a matrix porosity of 12.3%, and a matrix permeability of approximately [missing information]. The reservoir temperature was 75℃, and the formation water salinity was 85 g / L. Pore size distribution parameters obtained by combined cryogenic liquid nitrogen adsorption and mercury intrusion porosimetry experiments covered nanopores and micropores ranging from 3 to 800 nm, with fracture aperture distribution concentrated in the range of 50–300 μm. The Klinkenberg method was used to measure the pore pressure gradient. ~ The intrinsic permeability was obtained by fitting within the MPa / m range, and the calculated Knudsen diffusion coefficient of the nanopores was... The Fick diffusion coefficient of the micron-sized pore is The main diagonal components of the multi-scale effective transmission coefficient matrix after asymptotic expansion and assembly are shown in Table 1.

[0092] Table 1. Main diagonal components of the multi-scale effective transmission coefficient matrix

[0093]

[0094] In the discrete fracture network random generation stage, core fracture statistics and well logging imaging interpretation are used as joint constraints. Fracture trace length follows a log-normal distribution (mean 0.35 m, standard deviation 0.12 m), fracture aperture follows a gamma distribution (mean 150 μm, standard deviation 45 μm), and fracture dip and dip angle are constrained by microseismic location data to a dominant azimuth of 35° ± 15° north of east and a dip angle of 72° ± 8°. After 280 independent Monte Carlo sampling iterations, the sample statistical moments converge, generating an effective ensemble sample set. After 18 iterations of ensemble Kalman filter data assimilation, the posterior standard deviation of fracture parameters is reduced to aperture ± 12 μm and dip ± 6°. The principal components of the output equivalent permeability tensor are... direction μm², direction μm², the off-diagonal component reflects fracture anisotropy. ESS-100 bio-enzyme injection concentration was 500 mg / L, and ELS bio-enzyme injection concentration was 300 mg / L. Michaelis-Menten parameters for each component were determined by isothermal and isobaric core enzymatic reaction experiments at 75℃ and a mineralization of 85 g / L. Ethanol dehydrogenase component... for xylanase component for protease components for lipase components for Each component As shown in Table 2.

[0095] Table 2. Michaelis-Menten kinetic parameters of each enzyme component

[0096]

[0097] Adaptive multigrid enzyme kinetic operator splitting acceleration algorithm in 3D mesh ( The diffusion substep operates on a 5-layer algebraic multigrid, and the preset coupling residual threshold at the crack-matrix interface is set to [value missing]. The initial time step size of the enzyme-catalyzed reaction substep is set to After the first complete operator splitting cycle, the coupling residual at the crack matrix interface is as high as s. If the threshold is exceeded, the number of local mesh layers is increased from 5 to 6; after 12 complete iterations, the global coupling residuals are reduced to The following outputs the numerical solutions of the concentration fields at each time step of the grid nodes, provided that the preset thresholds are met. Figure 3 As shown, the concentration field distribution at different injection times indicates that the migration rate of the bioenzyme in the crack channel is significantly faster than that in the matrix pores. A clear concentration gradient transition zone is present at the matrix crack interface, which is consistent with the theoretical expectation of dual-pore media.

[0098] The enzyme concentration spatiotemporal evolution prediction model uses 300 sets of parameter combinations as training data, divided into training, validation, and test sets in an 8:1:1 ratio, with an initial learning rate of Adam optimizer. After 300 training iterations, the sum of the physical information loss on the validation set and the mean square error of the data fitting converged. In this test scenario, the concentration field adaptive control function value... If the value is greater than 0.8 in the first 5 time steps, it triggers an increase in the physical information loss weight coefficient and self-correction backtracking; starting from the 6th time step... Once the stability value drops below 0.4, the model enters a stable prediction phase. The weighting coefficients for physical information loss decrease, the mesh in uniform regions becomes coarser, and computational efficiency is significantly improved. In the final output of the three-dimensional concentration field time-series prediction results, the ESS-100 bioenzyme alcohol dehydrogenase component is mainly concentrated in the matrix nanopore region, while the ELS bioenzyme protease component is mainly distributed in the crack channels. The spatial distribution of both is highly consistent with the physical expectations of their respective dominant mass transfer mechanisms.

[0099] Compared to traditional single-scale deterministic fracture diffusion models, this invention uses multi-scale homogenization theory to mathematically connect the mass transfer coefficients of different pore scales, eliminating the mass transfer error caused by scale confusion when the Knudsen effect and fracture convection effect coexist in single-scale models. The integrated Kalman filter data assimilation method replaces analytical covariance with integrated statistics, eliminating the ill-conditioned nature of the Jacobian matrix in high-dimensional fracture parameter inversion from a mathematical structure perspective, enabling the equivalent permeability tensor estimation to continuously converge to the true state of the reservoir. The adaptive multigrid enzyme kinetic operator splitting acceleration algorithm removes rigid coupling constraints through Strang symmetric splitting and concentrates floating-point operations on the interface regions that require the most physical precision with an adaptive grid hierarchy strategy, ensuring global error tolerance while avoiding redundant calculations in uniform regions. The multi-scale spatiotemporal decomposition convolution-Transformer hybrid network embeds diffusion equation constraints into the training objective, enabling the artificial intelligence prediction model to still produce physically consistent concentration field time series results under sparse label conditions. Overall, it achieves multi-scale, multi-component, physically consistent closed-loop modeling capabilities that traditional numerical methods cannot achieve.

[0100] It should be noted that the variables involved in this invention are explained in detail in Tables 3 and 4.

[0101] Table 3. Variable Explanation Table (Part 1)

[0102]

[0103] Table 4. Variable Explanation Table (Part Two)

[0104]

[0105] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for establishing a bioenzyme diffusion model of tight reservoir matrix fractures, characterized in that, Includes the following steps: Core samples from tight reservoirs were collected, and pore-scale distribution parameters from nanopores to millimeter-scale fractures were obtained through nitrogen adsorption experiments and mercury intrusion porosimetry experiments. Based on the multi-scale homogenization theory, the Knudsen diffusion coefficient of nanopores, the Fick diffusion coefficient of micropores, and the fracture convection mass transfer coefficient were asymptotically expanded and uploaded to the Darcy scale to construct a multi-scale effective transport coefficient matrix of hierarchical transfer. The discrete fracture network random generation technique is used to construct integrated samples of fracture geometric topology for pore-scale distribution parameters. Dynamic pressure observation data is injected into the Monte Carlo forward propagation process through integrated Kalman filter data assimilation method to constrain the posterior probability distribution of fracture aperture, fracture orientation and fracture connectivity, and output the equivalent permeability tensor and matrix fracture mass exchange function. Using the equivalent permeability tensor and matrix fracture mass exchange function as input, the matrix fracture coupled diffusion reaction equation is decomposed into diffusion substeps and enzyme-catalyzed reaction substeps. An adaptive multigrid enzyme kinetics operator splitting acceleration algorithm is used to solve the problem iteratively on a multi-grid. After each complete operator splitting cycle, the coupling residual is calculated and the number of local grid layers and the number of smoothing iterations are dynamically adjusted. The numerical solution of the grid node concentration field at each time step is output. The multi-scale effective transport coefficient matrix, the equivalent permeability tensor and the matrix fracture mass exchange function, and the numerical solution of the concentration field of the grid nodes are input into the enzyme concentration spatiotemporal evolution prediction model. The enzyme concentration spatiotemporal evolution prediction model outputs the three-dimensional concentration field time series prediction results of ESS-100 bioenzyme and ELS bioenzyme in the matrix fracture dual-porosity medium. Using the coupling residual, the posterior standard deviation of the crack parameters, and the prediction residual of the three-dimensional concentration field time-series prediction results as inputs, the value of the concentration field adaptive control function is calculated. The physical information loss weight coefficient of the enzyme concentration spatiotemporal evolution prediction model is dynamically adjusted according to the interval to which the value of the concentration field adaptive control function belongs, and the adaptive multigrid enzyme kinetic operator splitting acceleration algorithm is re-driven to iterate until the global error tolerance is met. Based on the convergent three-dimensional concentration field time-series prediction results that meet the global error tolerance, and combined with the statistical results of the posterior probability distribution of crack parameters, the spatial distribution of each component of ESS-100 bioenzyme and ELS bioenzyme in the matrix and cracks is evaluated, and the final matrix crack bioenzyme diffusion model parameter set is output.

2. The method for establishing a biological enzyme diffusion model of tight reservoir matrix fractures according to claim 1, characterized in that, The multi-scale homogenization theory specifically uses asymptotic expansion as a mathematical tool to eliminate the microscopic heterogeneity in the pore-scale control equations through periodic cell integration, obtaining macroscopic effective parameters. Using the pore-scale distribution parameters as input, diffusion control equations are established for nanopores and micropores on the pore-scale cells, and periodic boundary conditions are applied. The Knudsen diffusion coefficient for nanopores and the Fick diffusion coefficient for micropores are obtained through numerical solution. After being transferred to the Darcy-scale equations through first-order asymptotic expansion, they together with the fracture convection mass transfer coefficients to form a multi-scale effective transport coefficient matrix.

3. The method for establishing a biological enzyme diffusion model of tight reservoir matrix fractures according to claim 2, characterized in that, The Knudsen diffusion coefficient of the nanopores was obtained by preparing core samples with a diameter of 25 mm and a length of 50 mm, conducting steady-state seepage experiments at 25 °C using nitrogen as the medium, and covering the pore pressure gradient. ~ At least 8 sets of flow rate data points were collected within the MPa / m range. The linear relationship between the inverse of pressure and apparent permeability was fitted using the Klinkenberg method. The intrinsic permeability was obtained by extrapolation to infinite pressure. The Knudsen diffusion coefficient of the nanopores was then calculated by combining the difference between the intrinsic permeability and the apparent permeability with the pore size distribution parameters.

4. The method for establishing a biological enzyme diffusion model of tight reservoir matrix fractures according to claim 3, characterized in that, The discrete fracture network random generation technique refers to using the Monte Carlo method to independently sample according to the multi-parameter probability distribution function of fracture trace length, fracture aperture, fracture dip direction, and fracture inclination angle, randomly generating a geometric set of fractures that satisfy statistical laws within the three-dimensional reservoir volume, and eliminating isolated fracture units through graph theory connectivity analysis, retaining the fracture network that constitutes the seepage channel; the multi-parameter probability distribution of fracture trace length, fracture aperture, fracture dip direction, and fracture inclination angle is jointly constrained by core fracture statistics, well logging imaging interpretation, and microseismic location data, and the parameter distribution is determined based on the convergence of sample statistical moments after at least 200 independent Monte Carlo samplings.

5. The method for establishing a biological enzyme diffusion model of tight reservoir matrix fractures according to claim 4, characterized in that, The integrated Kalman filter data assimilation method refers to approximating the error covariance in Bayesian updates with the integrated covariance matrix of the fracture parameter sample set. The residual between dynamic pressure observation data and forward model predictions is mapped to the fracture aperture, fracture direction and fracture connectivity parameter space through observation operators. Linear updates are performed synchronously for each integrated member, and the iteration continues until the statistical distribution of the fracture parameter sample set converges. The method outputs the posterior probability distribution of fracture parameters, the equivalent permeability tensor and the matrix fracture quality exchange function.

6. The method for establishing a biological enzyme diffusion model of tight reservoir matrix fractures according to claim 5, characterized in that, The adaptive multigrid enzyme kinetics operator splitting acceleration algorithm specifically decomposes the matrix crack coupled diffusion reaction equation into diffusion sub-steps and enzyme-catalyzed reaction sub-steps in each time step according to the Strang symmetric splitting strategy. The diffusion sub-step constructs a 5-layer grid hierarchy using an algebraic multigrid framework. In the finest layer, Gauss-Seidel iteration is used to eliminate high-frequency error components, and in the coarsest layer, a direct solver is used to correct low-frequency error components. The residuals and corrections are transferred between layers using standard constraint operators and extension operators. The enzyme-catalyzed reaction sub-step uses a Krylov subspace exponential time integrator to process the rigid nonlinear source terms of the Michaelis-Menten equations for each component of the ESS-100 and ELS bioenzymes.

7. The method for establishing a biological enzyme diffusion model of tight reservoir matrix fractures according to claim 6, characterized in that, ESS-100 bioenzyme is a complex system composed of protein, alcohol dehydrogenase and xylanase, which migrates in dense reservoir nanopores mainly through nanopore Knudsen diffusion. ELS bioenzymes are a complex system composed of proteases and lipases. The protease and lipase components migrate with the convective mass transfer in the crack channels and undergo enzymatic reactions on the crack walls. The alcohol dehydrogenase and xylanase components in ESS-100 bioenzymes, as well as the protease and lipase components in ELS bioenzymes, are parameterized by their respective maximum reaction rates and Michaelis-Menten constants in the Michaelis-Menten equations of the enzymatic reaction substeps.

8. The method for establishing a biological enzyme diffusion model of tight reservoir matrix fractures according to claim 7, characterized in that, The enzyme concentration spatiotemporal evolution prediction model is a multi-scale spatiotemporal decomposition convolution-Transformer hybrid network. The model input is a time-series snapshot sequence on a three-dimensional reservoir grid, with the multi-scale effective transmission coefficient matrix, equivalent permeability tensor, and numerical solutions of the concentration field of the grid nodes as channel features. The spatial dimension is designed with three levels of parallel convolution branches. The output feature maps of the three parallel convolution branches are uniformly upsampled to an intermediate resolution by a three-dimensional cross-scale feature pyramid fusion module and then stitched along the channel dimension. The stitched features are input along the time axis into a causal Transformer encoder. The encoder output is restored to the three-dimensional concentration field time-series prediction result by a fully connected decoding layer.

9. The method for establishing a biological enzyme diffusion model of tight reservoir matrix fractures according to claim 8, characterized in that, The specific steps for establishing the training dataset for the enzyme concentration spatiotemporal evolution prediction model include: using an adaptive multigrid enzyme kinetic operator splitting acceleration algorithm to train the dataset for a coverage porosity of 5%–25% and a crack density of 0.1–2.

0. 300 parameter combinations of ESS-100 and ELS bioenzymes with concentrations ranging from 10 to 5000 mg / L were run until convergence. A complete three-dimensional concentration field time-series snapshot of each combination was saved as a labeled sample. The training, validation, and test sets were divided in an 8:1:1 ratio. The training steps specifically included: using the Adam optimizer at an initial learning rate... The training process is repeated for 300 epochs. If the validation set loss does not decrease after every 50 epochs, the learning rate is reduced to 0.5 times the current value.

10. The method for establishing a biological enzyme diffusion model of tight reservoir matrix fractures according to claim 9, characterized in that, The physical information loss weight coefficient is dynamically adjusted according to the interval to which the concentration field adaptive control function value belongs. Specifically: when When the physical information loss weight coefficient is increased to 2.0 times its current value, a self-correcting iterative jump mechanism is triggered to recalculate the current time step; when At that time, the physical information loss weighting coefficient remains unchanged, and the number of local mesh layers at the crack matrix interface is increased by 1 layer only; when At that time, the physical information loss weight coefficient is reduced to 0.5 times the current value and the mesh is coarsened for uniform regions.

Images