A method for identifying response mode of water invasion and evaluating dynamic reserves of a fracture-type gas reservoir with edge and bottom water

Abstract

This invention discloses a method for identifying water intrusion response patterns and evaluating dynamic reserves in fractured edge-bottom water gas reservoirs, belonging to the field of oil and gas field development and reservoir engineering technology. It addresses the difficulties in distinguishing the coupling effects of water intrusion energy replenishment and gas containment, and the challenges in quantitatively evaluating dynamic water-sealed gas volume. The method acquires pressure, gas production, water production, and pressure-volume-temperature data to construct a dimensionless relative pseudo-pressure-production response curve. The response deviation is calculated using a closed gas reservoir baseline. Through historical data fitting with physical constraints, the net water intrusion energy replenishment coefficient, pressure reduction coefficient, and dynamic water-sealed gas formation parameters are obtained. The energy replenishment intensity, containment intensity, dynamic water-sealed gas volume, remaining movable gas volume, and mobility-restricted proportion are calculated to identify energy replenishment-dominated, containment-dominated, or energy replenishment-containment coexisting modes. This method can be used for dynamic diagnosis, reserve evaluation, drainage gas production, pressure-controlled production, and zoning adjustments in water-bearing gas reservoirs, improving the accuracy of pattern identification and reserve evaluation.

Description

Technical Field

[0001] This invention relates to the field of oil and gas field development and gas reservoir engineering technology, specifically to a method for identifying water intrusion response patterns and evaluating dynamic reserves in fractured edge-bottom water gas reservoirs. Background Technology

[0002] As the development of complex gas reservoirs, such as deep, ultra-deep, fractured carbonate rocks, and tight sandstones, deepens, water-bearing gas reservoirs, controlled by natural fractures, multi-scale heterogeneity, and edge / bottom water, often exhibit problems such as early water emergence, rapid water intrusion, steep decline, and low recovery rates in the later stages of development. In fractured edge / bottom water gas reservoirs, the fracture system is not only an advantageous channel for efficient natural gas seepage and production release, but also the main channel for edge or bottom water to preferentially intrude into the reservoir interior.

[0003] After rapid intrusion of edge-bottom water along highly permeable fractures, the increase in fracture water saturation, decrease in relative gas permeability, and reduction in local gas conductivity often occur simultaneously. At this time, a large amount of residual gas may still exist in the matrix, but its effective flow path with the wellbore is weakened or cut off by the water body. Some of the residual gas changes from a mobile state to a dynamic water-sealed state, resulting in a decrease in the utilization of the residual gas. Therefore, the key to dynamic analysis of fractured edge-bottom water gas reservoir development is not only to determine whether water intrusion is active, but also to identify the relative contributions of water intrusion energy replenishment and gas containment to the pressure-production response.

[0004] Existing water intrusion characteristic curve methods, water-gas ratio diagnostic methods, and development stage identification methods can reflect the intensity of water intrusion to some extent. However, most methods mainly focus on judging the level of water intrusion activity and cannot quantitatively distinguish the energy replenishment effect of net retained water on formation pressure and the sealing effect caused by the reduction of fracture conductivity. Existing mass balance methods can invert the cumulative water intrusion volume or dynamic water-sealed gas volume, but their focus is usually on volume calculation and it is difficult to directly generate energy replenishment-sealing response charts for rapid on-site identification.

[0005] Therefore, it is necessary to establish a method for dynamic interpretation of fractured edge-bottom water gas reservoir development, which unifies dimensionless relative pressure curves, closed gas reservoir baselines, response deviations, energy replenishment intensity indicators, and sealing intensity indicators to form identification criteria for energy replenishment-dominant, sealing-dominant, energy replenishment-sealing coexistence, and low-response boundary states, and further serve dynamic reserve evaluation and development adjustment. Summary of the Invention

[0006] The purpose of this invention is to provide a method for water intrusion response pattern identification and dynamic reserve evaluation in fractured edge-bottom water-gas reservoirs, which solves the problems that existing water intrusion diagnosis methods have difficulty distinguishing between water intrusion energy replenishment effect and gas containment effect, and that existing dynamic reserve evaluation methods lack response pattern identification charts and quantitative closed-loop of remaining movable gas volume.

[0007] This invention takes water intrusion replenishment-gas containment response mode identification as its core, and combines it with dynamic water-sealed gas volume, remaining movable gas volume and remaining gas mobility restriction ratio calculation to achieve the unity of response mode identification and dynamic reserve evaluation.

[0008] To achieve the above objectives, the present invention adopts the following technical solution; A method for water intrusion response pattern identification and dynamic reserve evaluation in fractured edge-bottom water gas reservoirs includes the following steps: S1: Obtain dynamic data on pressure, gas production, water production, fluid properties, and original geological reserves of the target gas reservoir. Calculate the dimensionless relative pseudo-pressure and recovery degree using the following formula, and establish... φ−R Develop response curves; S2: Using the closed gas reservoir baseline response as a reference, establish the closed gas reservoir baseline and calculate the target gas reservoir. φ−R The response curve deviates from the baseline, and the competitive response characteristics of water intrusion energy and gas containment are extracted based on the sign, amplitude, slope of change and zero-crossing position of the response deviation. S3: Based on target gas reservoir φ−R Response curves and response deviation curves, using the dimensionless quantity of purified water. I w,j The state variables to be determined are obtained by jointly inverting the net water intrusion energy coefficient through history fitting with physical constraints. c Pressure reduction factor or Dynamic water seal gas formation parameters a , b , d and the dimensionless quantity of purified water at each production moment. I w,j And calculate the energy replenishment intensity index. K e and enclosure strength index K s ; S4: Based on the energy replenishment intensity index K e The horizontal axis represents the sealing strength index. K s Using the vertical axis as the ordinate, a water intrusion energy replenishment-gas containment response intensity chart is constructed by combining at least two of the following: the positive and negative proportions of the response deviation, the maximum amplitude, the slope of change, and the zero-crossing characteristics. Based on the weak water intrusion-weak water containment verification samples, the total threshold for low response, the deviation amplitude threshold, the energy replenishment intensity threshold, and the containment intensity threshold are determined. S5: Calculated from the target gas reservoir K e , K sThe response deviation curve characteristics are matched with the response intensity chart to identify it as a type dominated by energy replenishment, a type dominated by isolation, or a type where energy replenishment and isolation coexist, and to identify whether it is in a low response boundary state. S6: Calculate the dynamic water seal gas volume, remaining movable gas volume, and remaining gas mobility restriction ratio at each production moment based on the parameter inversion results, and evaluate the dynamic reserve utilization status in conjunction with the model identification results to determine the development adjustment direction.

[0009] The pressure dynamic data includes average formation pressure, the gas production dynamic data includes cumulative gas production, the water production dynamic data includes cumulative water production, and the fluid property data includes natural gas deviation coefficient, initial natural gas volume coefficient, and formation water volume coefficient.

[0010] In this invention, the dimensionless relative pseudo-pressure, recovery degree, closed gas reservoir baseline, and response deviation are determined by the following formulas: ; ; ; ; In the formula, φ j For the first j The dimensionless relative pressure at each production moment; p j For the first j The average formation pressure at each production moment, in MPa; p i This represents the original formation pressure, expressed in MPa. Z j and Z i The first j Natural gas deviation coefficient at each production moment and initial state; R j To determine the degree of extraction; G p,j For the first j The cumulative gas production at each production moment, in cubic meters per second. 3 ; G Original geological reserves, in cubic meters (m³). 3 ; f 0, j The baseline response for a closed gas reservoir; Δ f j For the target gas reservoir φ−R The deviation of the response curve from the baseline of the closed gas reservoir; When Δ f j When Δ > 0, it indicates that the target gas reservoir pressure is maintained at a higher level than the baseline response of a closed gas reservoir, and the energy replenishment effect of the net retained water body is dominant; when Δ f j When Δ < 0, it indicates that the attenuation of gas phase conductivity and the reduction effect of gas sealing on pressure response are enhanced; when Δ < 0 between adjacent production times f j When a positive-to-negative transition occurs, it indicates that there is a phased competition between water infiltration energy and gas containment energy.

[0011] To ensure that the response chart identification is consistent with the mechanism of fractured edge-bottom water-gas reservoirs, this invention uses the dimensionless quantity of purified water intrusion. I w,j As state variables, they are solved sequentially according to production time, and the dynamic water seal gas ratio is adopted. S ws,j Characterizing the containment response; the quantities and the dimensionless relative pseudo-pressure calculation relationships are as follows: ; ; ; In the formula, I w,j For the first j The dimensionless quantity of clean water intrusion at each production moment; W e,j Cumulative water inundation, in meters (m). 3 ; W p,j Cumulative water production, in cubic meters (m³). 3 ; B w This is the formation water volume coefficient; B gi This is the initial natural gas volume factor; S ws,j For the first j The dynamic water seal gas volume at each production moment accounts for a portion of the original geological reserves. G The proportion; a Basic sealing factor; b To accelerate the sealing coefficient; d This is the marginal buffer coefficient; c The water replenishment energy coefficient; or This is the pressure reduction factor. In equation (5), the parameter inversion is first solved... I w,j Then convert W e,j .

[0012] In this invention, f j , R j , f 0,j Δ f j, I w,j , S ws,j , K e , K s and K tot All are dimensionless quantities. a、b、d、c and or All are dimensionless parameters; The objective function and parameter constraints for parameter inversion are as follows: ; ; ; In the formula, M This represents the number of historical fitted data points. f j obs For the first j The dimensionless relative pressure at each production moment is calculated from measured pressure data; f j cal For the first j The dimensionless relative pseudo-pressure calculated from the pressure response relationship at each production moment; l For constraint weights; P( i ) represents the physical constraint penalty item.

[0013] Parameter inversion employs a multi-parameter history fitting method with physical constraints. For any candidate parameter combination... θ= { a, b,d,c,h The pressure response closed-loop equation is solved according to the production time sequence. I w,j ;P( i ) by − I w,j , I w,j−1 - I w,j , d −(1−c I w,j ),−(1+d I w,j ),− S ws,j and S ws,j −(1− R j The sum of the squares of the positive parts of each quantity is constituted, where the positive part is the larger value of the quantity and zero, and the dimensionless quantity of the water intrusion before the first time step is set to 0.

[0014] When two mathematical roots exist, those that do not satisfy the condition are eliminated. I w,j ≥0、 I w,j ≥ I w,j−1 、1− cI w,j >δ、and 0≤ S ws,j ≤1− R j The root; when multiple physical solutions exist at the same time, choose the one that makes |I w,j −I w,j−1 | The smallest continuous root. d Since it is a dimensionless positive number, we can take 10. −8 Up to 10 −6 In this embodiment, 10 is used. −8 .

[0015] Parameter inversion can first use genetic algorithms, particle swarm optimization, differential evolution algorithms, or other global optimization algorithms to obtain candidate solutions, and then use nonlinear least squares methods, trust region methods, or sequential quadratic programming methods for local refinement. The inversion terminates when the relative change of the objective function between two consecutive iterations is lower than a preset convergence threshold, or when a preset maximum number of iterations is reached. In this embodiment, the convergence threshold is set to 10. −6 The maximum number of iterations is set to 2000.

[0016] This invention defines the diagnostic interval as a continuous data interval from the first valid production dynamic data point to the target evaluation time. Let the number of valid data points within the diagnostic interval be... N The average value of the state variables, the response intensity index, the total response intensity, and the maximum response deviation amplitude are determined by the following formulas: ; ; ; The low response total threshold, deviation amplitude threshold, energy replenishment intensity threshold, and sealing strength threshold of the weak water intrusion-weak water seal verification sample are defined as follows: ; ; ; ; In the formula, the superscript WW represents the weak water intrusion-weak water seal verification sample; N pos and Nneg Each of the diagnostic intervals satisfies Δ f j >0 and Δ f j The number of data points less than 0.

[0017] when K tot ≤ K tot,th ,and A Δ f ≤ A Δ φ,th When Δ is detected, it is identified as a low-response boundary state; after excluding low-response boundary states, when Δ f A positive-to-negative transition occurs within the diagnostic interval, or K e > K e,th and K s > K s,th When the energy replenishment-blocking coexistence type is not met, it is identified as a type where energy replenishment and blocking coexistence coexistence occurs; when no positive-to-negative conversion occurs and the identification conditions for energy replenishment-blocking coexistence are not satisfied, when... K e > K s and N pos > N neg When, it is identified as a supplementary energy-dominant type, when K s > K e and N neg > N pos When the condition is met, it is identified as a barrier-dominant type; otherwise, it is identified as a combination of energy replenishment and barrier.

[0018] The response intensity chart can be constructed by selecting at least two of the following: the positive and negative percentage of the response deviation, the maximum amplitude, the slope of change, and the zero-crossing feature; preferably, all four features are used for comprehensive identification.

[0019] The dynamic reserve evaluation described in this invention further calculates the dynamic water-sealed gas volume, the remaining movable gas volume, and the proportion of remaining gas mobility restricted: ; ; ; In the formula, G ws,j For the first j The amount of gas that is difficult to effectively participate in the flow at any given production moment due to water intrusion and blockage.G m,j The remaining gas volume still has effective utilization capacity. L m,j The proportion of the remaining gas that is confined by water intrusion; e To prevent the denominator from being zero, a very small positive number with the same dimensions as the gas volume is preferred. e =10 −12 G By constraining 0≤ S ws,j ≤1− R j ensure G m,j ≥0 and 0≤ L m,j ≤1.

[0020] The beneficial effects of this invention include: (1) The development response of fractured edge-bottom water-gas reservoirs is decomposed into two types of effects: water intrusion energy and gas containment. This can avoid misjudgment caused by judging the water intrusion type based solely on water production or water-gas ratio. (2) with φ−R Based on the response curve and the baseline of the closed gas reservoir, the response deviation is established, which can intuitively identify the dynamic characteristics of water intrusion energy replenishment, gas containment, and the phased competition between the two. (3) Adopt K e - K s The response intensity chart combines parameter sensitivity with actual water intrusion levels and dynamic water seal gas ratio, and uses a verification sample to determine the threshold, reducing the impact of relying solely on the original data. c and or Scale bias caused by absolute value partitioning pattern; (4) Combining response mode identification with dynamic water seal gas volume, remaining movable gas volume and the proportion of remaining gas mobility restriction can provide quantitative basis for judging the use of remaining gas, drainage gas extraction, pressure control production and zoning adjustment. (5) The data required by the method are mainly conventional pressure, gas production, water production and pressure-volume-temperature parameters, which are suitable for dynamic interpretation of fractured edge-bottom water gas reservoirs and mid-to-late stage reserve evaluation. Attached Figure Description

[0021] Figure 1 This is a schematic diagram of the overall process of the method of the present invention. Steps S1 to S6 correspond one-to-one with the method steps in the claims and specification. Figure 2 A schematic diagram of a dimensionless relative pressure-recovery response model and a baseline for a closed gas reservoir. Figure 3 For the response deviation Δ fDiagnostic diagram showing changes with the degree of extraction; Figure 4 Energy replenishment-isolation response intensity chart and verification diagram for four typical operating conditions; Figure 5 Different water replenishment energy coefficients c and pressure reduction factor or Dimensionless relative quasi-pressure response curves under the given conditions; Figure 6 For the field of the method of the present invention φ−R Historical fitting, Δ f - R Diagnostic and response intensity chart interpretation results. Detailed Implementation

[0022] The technical solution of the present invention will be further described below with reference to embodiments. These embodiments are used to illustrate the technical concept and application of the present invention and are not intended to limit the scope of protection of the present invention.

[0023] Example 1: Parameter Inversion and Response Intensity Chart Construction Pressure dynamics, gas production dynamics, water production dynamics, pressure-volume-temperature (PVT) data, and original geological reserves data of the target fractured gas reservoir with sidewater are selected. The dimensionless relative pseudo-pressure and recovery degree at each production moment are calculated according to formulas (1) and (2) to form the development response curve. The baseline of the closed gas reservoir is established according to formula (3), and the response deviation is calculated according to formula (4). In this embodiment, the search range for the parameters to be inverted is set as follows: 0≤ a ≤5、0≤ b ≤800、0≤ d ≤200、0≤ c ≤20、0≤ or ≤2. The above range is used for the inversion calculation in this embodiment and does not constitute a limitation on the range of parameter values ​​in this invention. During the inversion process, a global optimization algorithm is first used for multi-starting-point search, and then the top-ranked candidate solutions are locally refined; each set of parameters is solved step-by-step. I w,j The physical analysis was screened based on nonnegativity, monotonicity, positive denominator, dynamic water seal gas ratio, and nonnegativity of remaining movable gas. The diagnostic interval is taken from the first valid dynamic data point to the target evaluation time, and the average value of the state variables is calculated according to formulas (10) to (12). K e , K s , K tot and A Δ f Statistical data is simultaneously collected for each set of verification conditions. Npos , N neg And zero-crossing features, determined by the maximum value of the weak water intrusion-weak water seal check sample. K e,th =0.08 K s,th =0.06 K tot,th =0.14 and A Δφ,th =0.01, and the other three types of working conditions are used to check the energy replenishment-dominant area, the isolation-dominant area and the energy replenishment-isolation coexistence area respectively.

[0024] Example 2: Graphical verification of different water intrusion-water seal conditions A dual-pore, dual-permeability numerical model was used to provide pressure-production response samples. The model grid size was 40×40×8, with individual grid dimensions of 25m×25m×6m, corresponding to a model scale of 1000m×1000m×48m. The reservoir temperature was 115℃, the original formation pressure was 42.5MPa, the reservoir thickness was 48m, the matrix porosity was 0.09–0.18, and the matrix permeability was 0.01–0.1mD. In the model, the matrix system was mainly responsible for reservoir function, while the fracture system was mainly responsible for flow conduction and water intrusion channels, with high-permeability fracture zones set in the water intrusion direction. Based on the preset water intrusion input level and water seal parameter level, the verification samples are pre-set to four typical operating conditions: strong water intrusion-strong water seal, strong water intrusion-weak water seal, weak water intrusion-strong water seal, and weak water intrusion-weak water seal. The main parameter settings for the four operating conditions are shown in Table 1; Table 1 Parameter Settings for Four Types of Verification Conditions ; In Table 1, K a The aquifer permeability, K h and K v These represent the horizontal and vertical permeability of high-permeability fractures, respectively, and BHP. min For the lowest bottom hole flowing pressure, STG max This is the maximum daily gas production. This parameter setting makes the water intrusion input level and the water seal parameter level independent of each other, thereby avoiding the use of an undetermined output threshold to backdivide the verification sample; The verification results of the four types of operating conditions show that: in the strong water intrusion-weak water seal operating condition, the net water intrusion energy replenishment response is significant, the energy replenishment intensity index is relatively high, and the development response curve is mainly located above the closed gas reservoir baseline, which can be used as an energy replenishment-dominant sample; in the weak water intrusion-strong water seal operating condition, the dynamic water seal gas formation and pressure reduction response are strong, which can be used as a sealing-dominant sample; the strong water intrusion-strong water seal operating condition has both obvious energy replenishment and sealing responses, which can be used as an energy replenishment-sealing coexistence sample; the weak water intrusion-weak water seal operating condition has a small overall deviation and can be used to determine the low response boundary state; The relative errors of the cumulative water intrusion at the end of the four working conditions were +0.35%, +0.11%, -0.43%, and -0.60%, respectively, and the relative errors of the dynamic water seal gas volume at the end of the period were +0.49%, -7.43%, +9.27%, and +7.66%, respectively, indicating that the state quantities can be used for map verification and dynamic reserve evaluation. The diagrammatic identification relationships for different water intrusion-water seal conditions are as follows: (1) Strong water intrusion-strong water seal condition: The aquifer has high strength and the high-permeability cracks have a strong sealing response. K e and K s If all values ​​exceed the corresponding verification threshold, or if the response deviation changes from positive to negative, it is identified as a type where energy replenishment and isolation coexist. (2) Strong water intrusion-weak water seal condition: The aquifer has high strength and the crack sealing response is weak. K e > K s Furthermore, the response deviation is mainly positive, indicating that it is a compensation-dominated type; (3) Weak water intrusion-strong water seal condition: The total water intrusion scale is limited, and the high-permeability cracks provide strong sealing. K s > K e Furthermore, the response deviation was mainly negative, indicating a blockade-dominated type. (4) Weak water intrusion-weak water seal condition: Both the aquifer and the seal response are weak. K tot and A Δ f If none of them exceed the corresponding threshold, they are identified as low response boundary states.

[0025] Example 3: Field Application The syenite porphyry body of the Songkarsu Formation in the Dixi 18 well area of ​​the Kelamayi gas field was used as the field application object. The reservoir lithology in this area is dominated by syenite porphyry, and the reservoir space is mainly fracture-pore type, with an average porosity of 10.8% and an average permeability of 0.071 mD. Fractured reservoirs account for about 87.9%, and the gas-water interface elevation is about -3090 m. It belongs to a low-permeability fractured rimwater volcanic gas reservoir. Historical data was fitted using 152 on-site production dynamic data points, with a pressure range of 37.89 MPa to 28.26 MPa and a cumulative gas production range of 0.376 × 10⁻⁶. 8 m 3 Up to 3.983×10 8 m 3 The cumulative water production range is 0.04×10. 4 m 3 Up to 16.798×10 4 m 3 The results are obtained by inversion using formulas (1) to (9). a =1.4330, b =0.4141, c =3.0、 d =0、 or =2.0; the mean absolute percentage error (MAPE) of the dimensionless relative pseudo-pressure fitting is 1.15%, and the root mean square error (RMSE) is 0.0131; Calculated from each production moment on site I w,j and S ws,j Calculate according to formulas (10) to (12). K e =0.2475 and K s =0.2434, the response deviation changes from positive to negative when the extraction level is about 0.0238, therefore the area is identified as a type of coexistence of energy replenishment and isolation; Adopt a unified reserve caliber G =82.00×10 8 m 3 The cumulative gas production at the end of the period was 3.983 × 10⁻⁶. 8 m 3 The final dynamic water seal gas ratio is 0.224267; the dynamic water seal gas volume calculated according to formula (13) is 18.39 × 10 8 m 3 The remaining movable gas volume at the end of the period, calculated according to formula (14), is 59.627 × 10⁻⁶. 8 m 3 According to formula (15), the proportion of remaining gas with restricted mobility is 23.57%. In the final stage... I w,j =0.1500 condition, 1− cI w,j =0.5500>0 and S ws,j <1− R j The physical constraints are satisfied. Based on the identification results of this invention, the development and adjustment should take into account two aspects: on the one hand, control the ineffective water intrusion and prevent the edge and bottom water from further advancing along the high-permeability cracks; on the other hand, combine the dynamic changes of the remaining movable gas volume and the proportion of the remaining gas mobility restricted, and carry out drainage gas extraction, pressure control production or zonal adjustment for local isolation areas.

[0026] Without departing from the technical concept of this invention, the inversion algorithm, threshold determination method, data processing flow, and development and adjustment method can be replaced with equivalent ones.

Claims

1. A method for identifying water intrusion response patterns and dynamically evaluating reserves in fractured edge-bottom water-gas reservoirs, characterized in that, Includes the following steps: S1: Obtain dynamic data on pressure, gas production, water production, fluid properties, and original geological reserves of the target gas reservoir. Calculate the dimensionless relative pseudo-pressure and recovery degree using the following formula, and establish... φ−R Develop a response curve: ； ； In the formula, φ j For the first j Dimensionless relative pressure at each production moment; p j For the first j The average formation pressure at each production moment, in MPa; p i This represents the original formation pressure, expressed in MPa. Z j and Z i The first j Natural gas deviation coefficient at each production moment and initial state; R j To determine the degree of extraction; G p,j For the first j The cumulative gas production at each production moment, in cubic meters per second. 3 G represents the original geological reserves, in cubic meters (m³). 3 ; S2: Using the closed gas reservoir baseline response as a reference, establish the closed gas reservoir baseline and calculate the target gas reservoir according to the following formula. φ−R The deviation of the response curve from the baseline is used to extract the competitive response characteristics of water intrusion energy and gas containment based on the sign, magnitude, slope, and zero-crossing position of the deviation. ； ； S3: Based on target gas reservoir φ−R Response curves and response deviation curves, using the dimensionless quantity of purified water. I w,j The state variables to be determined are obtained by jointly inverting the net water intrusion energy coefficient through history fitting with physical constraints. c Pressure reduction factor η Dynamic water seal gas formation parameters a , b , d and the dimensionless quantity of purified water at each production moment. I w,j And calculate the energy replenishment intensity index. K e and enclosure strength index K s ; S4: Based on the energy replenishment intensity index K e The horizontal axis represents the sealing strength index. K s Using the vertical axis as the ordinate, a water intrusion energy replenishment-gas containment response intensity chart is constructed by combining at least two of the following: the positive and negative proportions of the response deviation, the maximum amplitude, the slope of change, and the zero-crossing characteristics. Based on the weak water intrusion-weak water containment verification samples, the total threshold for low response, the deviation amplitude threshold, the energy replenishment intensity threshold, and the containment intensity threshold are determined. S5: the calculated response intensity of the target gas reservoir is matched with the response intensity chart, and the type of the target gas reservoir is determined. K e , K s and the response deviation curve characteristics are matched with the response intensity chart, and it is identified as a dominant energy supplement type, a dominant sealing type or a coexisting type of energy supplement and sealing, and whether it is in a low response boundary state is identified. S6: Calculate the dynamic water seal gas volume, remaining movable gas volume, and remaining gas mobility restriction ratio at each production moment based on the parameter inversion results, and evaluate the dynamic reserve utilization status in conjunction with the model identification results to determine the development adjustment direction.

2. The method according to claim 1, characterized in that, In step S2, when Δ φ j When the value is greater than 0, the pressure response of the target gas reservoir is determined to be higher than the baseline response of the closed gas reservoir, indicating that the net retained water body has an energy replenishment contribution to the formation pressure. When Δ φ j When the pressure response of the target gas reservoir is less than 0, it is determined that the pressure response of the target gas reservoir is lower than the baseline response of the closed gas reservoir, indicating that the gas phase conductivity of the fracture is reduced and the residual gas sealing effect is enhanced; when the Δ at adjacent production times is less than 0, it indicates that the pressure response of the target gas reservoir is lower than that of the closed gas reservoir, indicating that the gas phase conductivity of the fracture is reduced and the residual gas sealing effect is enhanced. φ j A positive-to-negative conversion occurs, or φ−R When the response curve crosses the baseline of the closed gas reservoir, it is determined that there is a staged competition between water intrusion energy and gas containment.

3. The method according to claim 1, characterized in that, In step S3, the ratio of net retained water volume to the original formation natural gas volume is defined as the net water intrusion dimensionless quantity. I w,j The dimensionless quantity of purified water is then used as a state variable and solved sequentially according to the production time. j Cumulative water intrusion at each production moment W e,j Calculate using the following formula: ； In the formula, W p,j For the first j The cumulative water production at each production moment, in m³; B w This is the formation water volume coefficient; B gi This is the initial natural gas volume factor; G The original geological reserves; the net water intrusion dimensionless quantity is used to characterize the effective water intrusion intensity that has an effect on both pressure maintenance and gas containment.

4. The method according to claim 3, characterized in that, Dynamic water seal air ratio S ws,j The basic containment factor 'a' and the accelerated containment factor 'a' b and marginal buffer coefficient d The determination was made jointly, and the water replenishment energy coefficient was combined with the determination. c and pressure reduction factor η Establish a dimensionless relative quasi-pressure calculation formula, and calculate according to the following formulas: ； ； in, a Characterizes the initial containment capacity during the low-purity water intrusion stage. b Characterizing the degree of nonlinear enhancement in the formation of dynamic water-sealed gas in the middle and late stages, d Characterizes the marginal slowdown in the transformation of newly added clean water intrusion into dynamic water-sealed gas during the high water intrusion stage; c Characterizes the energy replenishment intensity of water intrusion on the mean pressure maintenance. η Characterizes the strength of the mapping from the dynamic water seal gas ratio to the average pressure reduction across the entire field.

5. The method according to claim 4, characterized in that, The parameter inversion aims to minimize the dimensionless relative pseudo-pressure history fitting error. The objective function and parameter non-negativity constraints are as follows: ； ； In the formula, M represents the number of historical fitted data points; φ j obs For the first j The dimensionless relative pressure at each production moment is calculated from measured pressure data; φ j cal The dimensionless relative quasi-pressure is obtained from the pressure response calculation formula at the j-th production moment; λ To constrain weights; P (θ) is the physical constraint penalty term. P (θ) is derived from − I w,j , I w,j−1 - I w,j , δ −(1− cI w,j ),−(1+d I w,j ),− S ws,j and S ws,j −(1− R j The sum of the squares of the positive parts of each quantity constitutes the sum of the squares of the positive parts of that quantity and zero, and the sum of the squares of the positive parts of each quantity before the first time step is given. I w ,0=0; Solve for each candidate parameter combination in the order of production time. I w,j When two mathematical roots exist, only those that simultaneously satisfy the condition are retained. I w,j ≥0、 I w,j ≥ I w,j−1 、1− cI w,j >δ、1+ dI w,j >0 and 0≤ S ws,j ≤1− R j The physical understanding, among which δ Preset to positive numbers; When multiple physical solutions exist, choose the one from the previous production time. I w,j−1 The consecutive roots with the smallest difference; The parameter inversion adopts a combination of global search and local optimization, and terminates when the relative change of the objective function is lower than a preset threshold or when the maximum number of iterations is reached.

6. The method according to claim 5, characterized in that, The diagnostic interval is a continuous data interval from the first valid production dynamic data point to the target evaluation time. Let the number of valid data points within the diagnostic interval be... N The average value of the state variables, the response intensity index, the total response intensity, and the maximum response deviation amplitude are calculated using the following formulas: ； ； ； ； ； ； ； in, N pos and N neg Each of the diagnostic intervals satisfies Δ φ j >0 and Δ φ j For data points less than 0, "WW" indicates a weak water intrusion-weak water seal verification sample; when K tot ≤ K tot,th ,and A Δ φ ≤ A Δ φ,th When Δ is detected, it is identified as a low-response boundary state; after excluding low-response boundary states, when Δ φ A positive-to-negative transition occurs within the diagnostic interval, or K e > K e,th and K s > K s,th When, it is identified as a coexistence of energy replenishment and containment; when no positive-to-negative conversion occurs and the identification conditions for coexistence of energy replenishment and containment are not met, when K e > K s and N pos > N neg When, it is identified as a supplementary energy-dominant type, when K s > K e and N neg > N pos When the condition is met, it is identified as a barrier-dominant type; otherwise, it is identified as a combination of energy replenishment and barrier.

7. The method according to claim 6, characterized in that, Based on preset water intrusion input levels and preset water seal parameter levels, the verification samples are pre-set as four typical operating conditions: strong water intrusion-strong water seal, strong water intrusion-weak water seal, weak water intrusion-strong water seal, and weak water intrusion-weak water seal, as the chart verification samples. Among them, the strong water intrusion-weak water seal condition is used to verify the energy replenishment-dominant area, the weak water intrusion-strong water seal condition is used to verify the isolation-dominant area, the strong water intrusion-strong water seal condition is used to verify the energy replenishment-isolation coexistence area, and the weak water intrusion-weak water seal condition is used to determine the low response total threshold and the deviation amplitude threshold. K e Threshold, K s Threshold and low response boundary state.

8. The method according to claim 4, characterized in that, The dynamic water seal gas volume G ws,j Remaining movable air volume G m,j and the proportion of residual gas mobility restricted L m,j Calculate according to the following formulas: ； ； ； In the formula, ε To prevent positive numbers with a denominator of zero, the constraint 0 ≤ S ws,j ≤1− R j ensure G m,j ≥0 and 0≤ L m,j ≤1; L m,j The larger the value, the higher the proportion of the remaining gas that is difficult to effectively participate in the flow due to water encroachment.

Images