An Integrated Multi-Scale Modeling Framework for Gas Entrainment Prediction in Coalbed Methane Production Systems

1 Introduction

As a clean and efficient unconventional energy resource, coalbed methane (CBM) is typically extracted via dewatering and pressure drawdown, a process designed to facilitate the desorption and subsequent recovery of methane from coal seams. Under ideal operating conditions, the pumping system is intended to lift only CBM produced water, while natural gas is recovered separately through the wellbore annulus. However, field monitoring across major global CBM basins—including the San Juan Basin (USA), the Bowen Basin (Australia), and the Qinshui Basin (China)—reveals that gas-liquid two-phase flow is prevalent within the production stream. This indicates significant gas entrainment in the produced water, deviating from the intended single-phase extraction [1,2,3].

The natural gas produced along with water mainly exists in two forms: first, dissolved methane, which is affected by reservoir temperature, pressure, and water chemical properties [4]. Second, dispersed bubbles, which are difficult to separate from water due to the negative pressure environment at the pump suction inlet [5]. Current investigations have shown that the natural gas entrainment in produced water from a single well can reach 0.5–3.2 m3/d. Based on the estimation of the global annual CBM production of 70 billion m3, the annual loss of natural gas resources along with produced water exceeds 3 billion m3 [6,7]. Moreover, most unrecovered natural gas is released directly into the atmosphere during produced water treatment. Because methane’s 100-year GWP is 28 times higher than that of CO2, its emission is a major source of greenhouse gas pollution [8,9]. A study by Omara et al. of over 1000 sites found that episodic events, particularly wellhead fluid discharge, are the leading cause of emissions in gas fields [10]. Similarly, studies at the Shengli Oilfield show associated gas emissions from produced water reaching 0.415 g/(m3·h) [11].

The current research on tools used for water drainage during coalbed methane (CBM) extraction emphasizes the development and application of various pump technologies, including rod pumps and rod-less alternatives. Traditional rod pumps remain a significant component in CBM drainage, with studies highlighting their operational mechanisms and technological advancements [12]. These pumps are widely used due to their reliability and effectiveness in lifting water and gas from coal seams, although challenges such as equipment durability and operational efficiency persist.

To elucidate the mechanisms governing gas ingress into produced water pipelines, we delineate the fluid transport trajectory from the coal seam to the wellbore (Fig. 1). Initially, the millennial-scale coexistence of water and methane under reservoir pressure establishes a state of saturated solubility. As extraction-induced pressure drawdown occurs, this dissolved gas exsolves into a free phase. During the recovery process, fluctuations in the wellbore liquid level or the positioning of production intervals below the pump intake—typically a reciprocating rod pump—can trigger the mechanical suction of gas. Furthermore, the structural integrity of wellbore tubulars is frequently compromised by recurrent mechanical wear and corrosion; such fractures not only facilitate internal fluid leakage but also permit the infiltration of extraneous gas from the surrounding formation. Consequently, characterizing the gas content in produced water presents a multi-faceted challenge. Current research is largely confined to regional assessments based on atmospheric gas escape, leaving a critical void in well-specific, high-resolution diagnostics. While high-fidelity Computational Fluid Dynamics (CFD) offers sophisticated insights into gas-liquid multiphase flow, the inherent complexity of 3D numerical simulations—compounded by dynamic boundary conditions and phase transitions—results in prohibitive computational latencies, often spanning several weeks for a single solution [13,14,15]. This bottleneck precludes real-time prognostic assessment of volatile downhole conditions and hinders the fundamental investigation of gas-bearing mechanisms. To bridge this gap, reduced-order modeling (ROM) and data-driven reconstruction techniques have emerged as transformative paradigms, offering a robust balance between physical rigor and computational efficiency in resolving complex multiphase dynamics [16,17].

In addition, for the recovery of natural gas from produced water, it is difficult to carry out separation and recovery directly at the wellhead due to the high cost. Actually, gas-liquid mixtures are usually transported from individual wells to processing stations through pipelines for centralized separation. During pipeline transportation, the flow regime of the gas-liquid mixture changes with the ups and downs of the pipeline route, resulting in unstable transportation flow rates, and fluctuations in the gas-liquid flow rates entering the centralized separators, even slug flow. These phenomena will have a significant impact on the design of the gas-liquid separator [18].

Dottore et al. [19,20] elucidated the pump performance dynamics in gas-bearing wells by formulating analytical expressions for pump chamber pressure fluctuations, volumetric efficiency, and the critical thresholds for gas lock initiation. From a mechanistic perspective, Das et al. [21,22] investigated gas–liquid flow transitions within vertical annular pipes, identifying gas velocity fluctuations at constant pressure as the primary driver of regime shifts. Furthering this, Gomez et al. [23] developed a unified steady-state mechanistic model to predict flow regimes, liquid holdup, and pressure gradients. Despite these advancements, a significant knowledge gap persists: existing frameworks seldom integrate the complex interplay between two-phase pipeline pulsations and the transient flow phenomena within gas–liquid separators.

In recent years, with the development of techniques such as order reduction and machine learning, scholars have begun to conduct large-scale model predictions based on well measured data and mathematical models for pipe pumps. Jamalbayov et al. [24] adopted a discrete modeling approach to adjust and decouple complex control equation systems, proposing a comprehensive model for the development of rod pump wells and reservoirs. Zuo et al. [25] based on a large amount of measured data on the electric power of rod pumps, used the MSISSA-BiLSTM learning model to estimate effective liquid production, forming an efficient prediction method. Therefore, techniques such as order reduction and machine learning have strong application potential in this type of problem.

A certain coalbed methane field mainly uses rod pumps to extract produced water. In recent years, the gas field has been plagued by the problem of gas content in produced water. To address the gas-entrainment challenges in CBM produced water, this study develops a multiscale integrated framework combining High-Fidelity CFD, Reduced-Order Modeling (ROM), and Ensemble Learning. We first conduct 3D two-phase flow simulations to characterize the transport mechanisms within single-well extraction systems. To overcome computational bottlenecks, these 3D physics are distilled into an efficient 1D numerical model via ROM, significantly accelerating individual analysis cycles. Leveraging this efficiency, a Stacking-based ensemble model is constructed to predict well-scale gas-liquid flow rates and gas content across the entire extraction block. Furthermore, we establish a network-scale two-phase flow model to analyze the impact of topographical undulations on transport characteristics. Finally, the separation performance of tubular gas-liquid separators is evaluated under realistic pulsating inflow conditions using 3D numerical analysis, providing a holistic understanding of the gas trajectory from the reservoir to surface processing.

Figure 1: Schematic diagram of gas source in produced water.

2 Research Methodology

2.1 Numerical Analysis Method

Rod pumps are commonly used for extracting produced water in the studied CBM extraction area. In the analysis of downhole produced water flow, it is necessary to consider the displacement and mechanical properties of the pump rod, piston, and valve ball, while investigating the impact of the suction effect generated by the operation of the pump and pipeline. In three-dimensional flow field CFD, the governing equations mainly include the mass, the momentum and the energy conservation equation: Mass conservation equation:∂ρ∂t+∇·ρV=S(1) Momentum conservation equation:∂ρ·V∂t+∇·ρVV=−∇P+∇μ∇V+∇VT+ρg+Tσ(2) Energy conservation equation:∂ρT∂t+∇·ρVT=∇·λ∇T(3) where: ρ is density, t is time, ∇ is the differential Hamiltonian operator, V is the velocity vector, P is pressure, S is the mass source term, μ is viscosity, T is temperature, g is gravitational acceleration, and T σ is the surface tension at the phase interface, λ is the equivalent thermal conductivity of the fluid medium. All physical quantities are mixed-phase parameters. The compressibility of the gas is using the ideal gas state equation.

Since the natural gas in the studied gas field has an extremely high methane content, the solubility of natural gas in produced water can be calculated using the following formula. Furthermore, the amount of dissolved/precipitated gas caused by pressure changes can be derived in the numerical model: x=fCH4gγx·Hx0(4) where: f CH 4 g is the fugacity of methane in the gas phase, calculated using the gas state equation. γ x is the activity coefficient of methane in water. and H x 0 is the reference-state Henry’s law constant.

The produced water in the studied gas field exhibits a salinity of 1760 mg/L, a measured viscosity of 1.02 mPa·s, and a surface tension of 76 mN/m.

Since the flow channel changes in CFD calculations are mainly dominated by the overall lifting and lowering of the plunger and valve ball. Therefore, it is assumed that only the fluid-structure interaction interfaces of the plunger and valve ball deformation during calculations, and the deformation mode is very similar to one-dimensional elastic deformation. Thus, based on Hooke’s law, grid layers near the wall are defined as the deformation zone, and the force balance relationship of each grid node in the deformation zone can be expressed by the following formula: Fi=∑jniKijΔXj−ΔXi(5) where: ni is the number of adjacent nodes of node i. K ij = 1 X i − X j is the spring stiffness coefficient between nodes i and j. Δ X i and Δ X j represent the displacements of grid nodes i and j relative to their initial positions.

Once the displacements of nodes on the fluid-structure interaction interface are obtained, the displacements of nodes in the flow channel are determined by the following formula: Xin′+1=Xin′+γΔXim(6) where: n′ is the time step. X i n ′ and X i n ′ + 1 are the positions of node i at the current time step and the next time step. γ is the relaxation factor, which is used to characterize the of the displacement of the node. Δ X i m is the displacement of node i at the current time step. and m is the number of iterations.

The above expressions for the natural gas dissolution/precipitation process and fluid-structure interaction process are implemented in the computational model through C language subroutines in CFD calculations, respectively applied to the mass conservation equation and the momentum conservation equation.

The downhole numerical analysis domain includes the area enclosed by the wellbore, plug, pump, and tubing. A three-dimensional geometric model of the pipeline and rod pump was established using 3D CAD methods, and the 3D analysis model generated by the filling method is shown in Fig. 2a. The analysis model of the gas-liquid separator is shown in Fig. 2b.

The flow field analysis model was generated using filling technology to ensure consistency with the solid model at the interface. Structured hexahedral (or quadrilateral) topological structures were used for meshing. The grid-independence of the numerical solution for the temperature was examined to ensure the grid independent solution. By varying the grid scale, 5 grid systems were generated. For the under well model, taking the conditions of pump setting depth of 817 m, pump length of 8.3 m, and pump submergence of 1 m. For the gas-liquid separator model, taking the conditions of an inlet liquid flow rate of 34.25 m3/h, and a gas flow rate of 27.29 m3/h. The results of the last two sets of grids are fairly close to each other (shown in Fig. 3), and the error was less than 1%. This indicates that the resolution of the finer grid was adequate to guarantee computational precision. The final grid number used for under well simulations was 5.54 million, and for separator simulations was 6.27 million.

The detailed model selection is shown in Table 1 [26,27,28,29]. In the selection of mathematical models, reference was made to existing research on gas-liquid two-phase flow in vertical pipelines [30]. To ensure sufficient computational accuracy, the selection of the difference scheme was determined by referring to the CFD calculation recommendations of AIAA [31]. In the context of vertical pipeline flow, the VOF model has been extensively employed to simulate and analyze gas-liquid interactions [32]. Therefore, the VOF model was employed in the calculation to simulate two-phase flow.

All boundary conditions for numerical calculations are derived from actual under well sensor measurement data. In the CFD calculations, no less than 3000 time iterations were done to reach a fully converged solution: Stability of all monitoring parameters (flow rate, force), mass balance convergence and maximal residuals about 10−5.

Figure 2: 3D numerical analysis model, (a) Under well, (b) Gas-liquid separator.

Figure 3: Grid independence assessment, (a) Under well, (b) Gas-liquid separator.

2.2 Reduced Order Analysis Method

Reduced-order model (ROM) analysis utilizes input-output data from a numerical analysis system to determine a simple model system through a type of complex model system, making the two equivalent in terms of solution capability [33,34]. In this study, the complex model is a downhole numerical analysis based on 3D transient numerical calculations.

In downhole numerical analysis, the displacement of the valve ball under solid constraints and two-phase countercurrent impact is the input of the system, and the gas-liquid flow rate at the pump outlet is the output of the system. To establish a one-dimensional transient numerical calculation model using this method, it is first necessary to obtain the hydraulic coupling dynamic characteristics of the valve ball. The purpose of model order reduction is to obtain this mechanical characteristic model. The system identification process involves first determining the design of the identification scheme and model type, then determining the parameters of the system based on characteristics of the system, and finally obtaining the model parameters, as shown in Fig. 4.

Figure 4: Reduced order analysis process diagram.

Parameter estimation of the model is a key step in the system identification process, which involves using iterative estimation methods to estimate model parameters to minimize the likelihood function. Based on the transient gas-liquid two-phase flow dynamic reduced-order model derived from system identification, the hydraulic coupling state equation of the valve ball can be expressed by the following formula. For a general fluid-structure interaction elastic system: x˙at=Aaxat+Bbξtfat=Ccxat+Ddξt(7) where: x a is the state vector. f a is the two-phase flow dynamic coefficient. ζ = u is the modal displacement. and A a , B b , C c , D d is the coefficient matrix to be identified.

In the solution process of 3D transient numerical calculations, the displacement of the valve ball can be taken as the input signal, and the corresponding two-phase flow force is obtained as the output signal. Then, system identification technology is used to perform parameter identification on the system. Therefore, for a given calculation condition, only one input excitation signal needs to be calculated to obtain the reduced-order identification model [35].

The study adopts the Proper Orthogonal Decomposition (POD) based on multivariate statistics for system parameter identification. Its basic process is to find a new set of bases in a large multivariate system, such that the dimension of the new system spanned by these new bases is significantly reduced, and the error between the new system and the original system is minimized. Specifically, if the solution obtained by solving the required partial differential equations is Y = y 1 , ⋯ , y p , where y i are column vectors, this solution matrix can be used as snapshots. A set of eigen orthogonal bases Ψ = φ 1 , ⋯ , φ p is derived from these snapshots, and the original system is projected onto this set of orthogonal bases, with the error between the projected system and the original system minimized.

The computational load of the hydraulic coupling response using the reduced-order method is far less than that of the direct 3D transient numerical calculation method. After multiple identifications, the identified hydraulic coupling dynamic model of the valve ball can be used to calculate the gas-liquid two-phase flow output rate of the pump.

3 Flow Field Characteristics in Wells

3.1 Flow Field Analysis

During the extraction process, the pump inlet may contain a certain proportion of gas because it is close to the coal seam or the liquid level outside the pump inlet is low, leading to gas being carried into the pump. This is the main cause of gas content in produced water. Taking the typical condition where gas overflows from the underlying reservoir at the pump inlet, a flow field analysis was conducted. The flow field distribution in the wellbore annulus below the pump inlet is shown in Fig. 5. During the upstroke of the pump, due to the pump’s suction effect, the gas entering the annulus through the perforations forms a regular gas column distribution and ultimately flows to the pump’s suction inlet. In the downstroke, however, the driving force for forced flow in the annulus stops, and the gas then presents an irregular random distribution. Nevertheless, it still moves upward overall under the drive of buoyancy. The flow field near the pump is shown in Fig. 6. Due to the suction effect of negative pressure, gas is entrained into the pump, causing gas carry-over issues.

Figure 5: Distribution of 3D flow field in the annulus.

Figure 6: Phase distribution under the condition of gas at the pump suction inlet.

The time-history curves of gas and liquid discharge flow rates after the pump stabilizes over multiple working cycles are extracted, as shown in Fig. 7a. Due to the gas content at the suction inlet, the pump chamber cannot continuously and stably inhale liquid, resulting in a significant reduction in liquid volume. From the transient gas-liquid flow rate curves, although the gas volume under standard conditions is much higher than the liquid volume, the liquid volume fraction in the fluid still dominates under downhole pressure.

As an important input for reduced-order analysis, the displacement time-history curve of the valve ball is shown in Fig. 7b. During the upstroke, as the piston moves upward, the traveling valve ball closes under the fluid pressure in the annular space, blocking the backflow of fluid in the pump barrel. Meanwhile, the negative pressure inside the pump barrel pushes opens the fixed valve ball, allowing fluid to enter the pump barrel. During the downstroke, the piston moves downward to squeeze the fluid in the pump barrel, and the increased pressure pushes open the traveling valve ball. The fluid then enters the inner channel. The pressure inside the pump barrel is higher than that in the annular space, so the fixed valve ball closes under the hydraulic pressure to prevent fluid backflow to the casing-tubing annulus. After each valve ball reaches the maximum or minimum displacement point, its displacement should theoretically remain unchanged for a short period. However, the effect of two-phase flow on the ball is far less stable than that of single-phase incompressible liquid, and slight oscillations of the ball are clearly observed during this process.

Using the reduced-order analysis method, after multiple iterations to obtain the hydraulic coupling dynamic model of the valve ball and applying it to the 1D model, the calculation results are very close to those of the 3D calculation. The total calculation time is reduced from nearly 200 h to 27 min.

Figure 7: Time history distribution of flow rate and valve ball displacement, (a) Gas liquid flow rate, (b) Valve ball displacement.

To evaluate the accuracy of the reduced-order model, comparative calculations were performed under nine different well conditions, and the proposed 1D reduced-order model was compared with the 3D model. The comparison results for gas and liquid flow rates are shown in Fig. 8, indicating that the two models yield very close results. The root mean squared error (RMSE), mean absolute error (MAE), and coefficient of determination (R2) were selected as the prediction error evaluation metrics. For the 1D reduced-order model, the RMSE, MAE, and R2 values for gas flow rate prediction are 5.55 L/min, 1.77 L/min, and 0.96, and for liquid flow rate prediction are 0.18 L/min, 0.21 L/min, and 0.95, respectively.

Figure 8: Comparison between 3D and 1D well flow rate calculation, (a) Gas flow rate, (b) Liquid flow rate.

3.2 Predictive Large Model

To ensure model parsimony while bolstering predictive performance, we developed an ensemble framework that integrates heterogeneous base learners to capture diverse data features. Specifically, Support Vector Regression (SVR) [36], Back-Propagation Neural Networks (BPNN) [37], and Random Forests (RF) [38]—representing kernel-based, neural-based, and tree-based architectures, respectively—were selected as base learners. These were coupled with an Adaptive Neuro-Fuzzy Inference System (ANFIS) as the meta-learner [39]. This stacking ensemble strategy [40] facilitates the construction of a robust proxy model for gas volume prediction, the architecture of which is illustrated in Fig. 9. The implementation follows a tripartite process:

• (1)Dataset Generation Stage

The sample dataset (X, Y) is divided into a training dataset ( X tr , Y tr ) containing numerical analysis data and a test dataset ( X te , Y te ) containing on-site measured data. X={xi,i=1,2,3…,n}(8) where: x i represents the input parameters of the i-th dataset sample, including 12 parameters: casing diameter, casing pressure, pump setting depth, dynamic liquid level depth, pump length, pump diameter, stroke length, stroke frequency, screen length, screen outer diameter, total gas production of the underlying reservoir, average depth of the underlying reservoir. Y=yi,i=1,2,3…,n(9) where: y i represents the response value corresponding to the j-th dataset sample, including the gas flow rate and liquid flow rate. n − k samples are selected as the training dataset:

According to the on-site conditions of the field, n − k = 317 samples are selected as the training dataset, k = 120 samples are selected as the test dataset:

• (2)Training Stage

Five-fold cross-validation is adopted to train each base learner (SVR, BPNN, RF). For each base learner, the original training set ( X tr , Y tr ) is first randomly and equally divided into five subsets (Training Sets 1–5). In each iteration, four subsets are used as the training set to train the base learner, and the remaining one is used as the validation set to evaluate the model’s prediction performance and obtain a prediction result Ptr−j. This imposes an implicit regularization on the base learners, which greatly reduces the risk of overfitting. This process is repeated five times to obtain the prediction results of the base learner for the entire original training set. For example, for the SVR learner:

The above training process is performed for each base learner to obtain a new training set ( P tr , Y tr ) constructed by the prediction results of each base learner and the corresponding original response values. This new training set is used as the training samples for the ANFIS meta-learner, where:

• (3)Testing Stage

The trained base learners are used to obtain the predicted values P te corresponding to the original test set. A new test set ( P te , Y te ) is constructed by these predicted values and the response values in the original test set, which is used to test the trained meta-learner, where:

Thus, the final prediction results F te = { f i , i = 1, 2, 3…, k} are obtained. Based on the above steps, a Stacking ensemble learning surrogate model for gas content prediction can be constructed. To improve the prediction accuracy of the surrogate model, the improved sparrow search algorithm (ISSA) is further adopted to optimize the parameters of each learner in the model.

Figure 9: Integrated learning framework.

Using the reduced order model described in the previous section, the learning samples can first be generated. A total of 317 sets of different well operating conditions are produced as training datasets. Each set of data contains 12 input parameters (casing diameter, casing pressure, pump setting depth, dynamic liquid level depth, pump length, pump diameter, stroke length, stroke frequency, screen length, screen outer diameter, total gas production of the underlying reservoir, average depth of the underlying reservoir) and 2 output parameters (gas flow rate, liquid flow rate).

On this basis, operational testing was conducted on 72 wells in the field to obtain gas and liquid flow rates at the wellhead as validation datasets. To demonstrate the accuracy of the developed prediction model, the proposed Stacking learning model was compared with individual base learners (SVR, BPNN, and RF) based on the gas flow rate obtained from validation datasets. The root mean squared error (RMSE), mean absolute error (MAE), and coefficient of determination (R2) were selected as the prediction error evaluation metrics. The comparison results of the four models are shown in Fig. 10.

All four models are capable of predicting gas flow rate with reasonable accuracy. Compared with the other three learners, the Stacking learning model achieves the highest prediction accuracy, with RMSE, MAE, and R2 values of 1.75 L/min, 0.82 L/min, and 0.98. This indicates that the model can complement the strengths of each base learner during the data learning process, demonstrating excellent predictive capability, and performing more robustly and accurately.

Figure 10: Comparison chart of prediction error.

4 Flow Field Characteristics of Transportation Separation System

4.1 Produced Water Pipe Network

The coalbed methane field under study is in a hilly terrain, with produced water pipelines laid along the surface showing obvious undulations. The entire field includes 101 wells distributed over an area of approximately 50 km2. A 1D pipeline network model was established as shown in Fig. 11, which marks the average annual water production of each production point. The pipe material is PE, and the calculation utilizes the average scaling thickness of 3.63 mm and surface roughness of 0.27 mm obtained from on-site investigations. The pipeline network is distributed in a dendritic pattern, converging at the central treatment station. The prediction model obtained in Section 3.2 was used to calculate the water and gas flow rate of each well in a typical year, with the calculation conditions selected as the maximum daily value, minimum daily value, and annual average value of the total regional water production in a year. The two-fluid model is employed to simulate 2 phase flow in pipelines, and the flow pattern state is determined using a region identification method based on flow pattern maps [41].

The pressure at the inlet of each production point’s pipeline is marked in Fig. 9. Since the outlet pressure of the pipeline network at the treatment station is close to atmospheric pressure, the pressure at the inlet of each production point’s pipeline represents the pressure drop along each pipeline route. The gas holdup entering the treatment station is shown in Fig. 12. From Fig. 12a, it can be found that the gas holdup fluctuates significantly, and there are moments when it is entirely gas or entirely liquid. In addition, the influence of the maximum or minimum water volume on the gas holdup is not obvious from a time-history perspective.

Fast Fourier Transform (FFT) transformation was performed on the gas holdup time-history data to obtain the pulsation distribution of gas holdup at different frequencies, as shown in Fig. 12b. At low frequencies, the pulsation energy under different water production conditions shows one or more peaks, while the pulsation energy decreases rapidly with the increase in frequency. The peak of the maximum pulsation frequency ranges from 0.008 to 0.007 Hz, corresponding to a pulsation period of 120–150 s. This indicates that the gas-liquid pulsation at the outlet of the pipeline network is dominated by low-frequency pulsation, which will have an adverse impact on the operation of the gas-water separator planned to be installed at the treatment station.

Figure 11: Numerical analysis model for extracted water pipe network.

Figure 12: Distribution of gas volume fraction at the end of the pipeline network, (a) Time History, (b) Frequency history.

4.2 Gas-Liquid Separator

This study adopts a tubular separator as the separation equipment for the centralized treatment of produced water. To verify whether the designed separator can operate properly—especially to well accommodate the incoming medium under significant gas-liquid flow rate pulsation. The annual maximum water production condition obtained in Section 4.1 was chosen as the calculation condition for the separator. The time-history gas and liquid flow rates calculated from the pipeline network were used as the inlet boundary conditions for the separator, and the internal flow field distribution under stable working conditions is shown in Fig. 13.

After entering through the separator inlet, the gas-liquid mixture is directed to 8 parallel branch pipes. As the flow channel space expands in the branch pipes, the flow velocity decreases rapidly. Reducing the flow velocity can effectively reduce the inertial force of the fluid, allowing the gas to be gradually dominated by buoyancy. Since the gas phase outlet is located at the highest point of the entire equipment, the gas can flow out smoothly from this outlet. Meanwhile, the liquid continues to flow downstream through branch pipes at the bottom. During operation of the separator, stable stratified flow can be maintained in the branch pipes, ensuring that the downstream liquid phase outlet is always below the liquid level and preventing gas from flowing in. From the perspective of flow field distribution, the separator can maintain a clear gas-liquid interface even under the conditions of maximum flow rate and strongest gas-liquid pulsation, and the interface level does not exceed the mid-plane of the lower branch pipes.

The liquid and gas holdup at the outlet of the separator are extracted as shown in Fig. 14. The liquid holdup at the gas outlet is slightly higher than the gas holdup at the liquid outlet, indicating that a small amount of liquid still fails to fully settle during the flow process. However, both the liquid content at the gas outlet and the gas content at the liquid outlet are maintained below 12‰, the capture efficiency of liquid droplets with a particle size larger than 150 μm can reach 99%, demonstrating that the separator has a good separation effect without obvious gas/liquid carry-over.

Figure 13: Flow field distribution of gas-liquid separator, (a) Streamline diagram, (b) Liquid phase distribution.

Figure 14: Distribution of phase content at the outlet of gas-liquid separator.

5 Conclusion

This study elucidates the complex dynamics governing gas entrainment in coalbed methane (CBM) produced water through a multi-scale framework that integrates high-fidelity 3D-to-1D reduced-order modeling (ROM) and advanced ensemble learning. By projecting complex 3D fluid dynamics onto a computationally efficient 1D manifold, we have developed a framework capable of rapid, cross-well characterization, providing a definitive depiction of gas entrainment processes that were previously computationally prohibitive to resolve.

Mechanistically, our analysis identifies the ‘suction effect’ of the rod pump as the primary driver for gas ingestion; specifically, intense localized flow gradients formed during the upstroke facilitate the entrainment of gaseous phases from the suction inlet into the pump chamber. To scale these insights for field-wide application, we established an optimized Stacking ensemble learning architecture—synthesizing SVR, BPNN, and Random Forest algorithms via an ANFIS meta-learner—which achieves remarkable predictive fidelity, with RMSE, MAE, and R2 values of 1.75 L/min, 0.82 L/min, and 0.98.

Furthermore, systemic investigation of the surface infrastructure reveals that while the gathering network is dominated by low-frequency pulsations, the peak of the maximum pulsation frequency ranges from 0.008 to 0.007 Hz. The downstream tubular separators possess sufficient inertial buffering capacity to maintain a stable gas-liquid interface, both the liquid content at the gas outlet and the gas content at the liquid outlet are maintained below 12‰, the capture efficiency of liquid droplets with a particle size larger than 150 μm can reach 99%. This comprehensive modeling suite not only resolves the mechanistic origins of gas-bearing produced water but also provides a robust tool for the real-time prognostic assessment and optimization of CBM extraction and transportation systems.