An Analytical Model Coupled with Orthogonal Experimental Design Is Used to Analyze the Main Controlling Factors of Multi-Layer Commingled Gas Reservoirs

: The majority of China’s multi-layer low permeability tight gas reservoirs are currently being extracted through the method of multi-layer co-production. However, due to the signiﬁcant disparity in physical properties and varying degrees of pressure depletion among the production layers, elucidating the primary factors inﬂuencing the productivity contribution of each gas layer remains challenging. A multi-factor analytical model is proposed for commingled gas wells with multiple layers. An unstable model is established for the production of commingled layers, and the problem of ﬂow distribution is addressed using the Duhamel convolution principle. The Laplace transform is subsequently employed to derive the solution in the Laplace domain, which can be inverted utilizing the Stehfest inversion algorithm to obtain a real-time domain solution. The inﬂuence of reservoir factors on the stratiﬁcation contribution rate has been comprehensively analyzed, encompassing permeability, porosity, initial pressure, drainage radius, and layer thickness. The orthogonal test design was employed to conduct range analysis and variance analysis separately, yielding the primary and secondary order as well as inﬂuence weight of the ﬁve factors. The ﬁndings demonstrate that, within this gas reservoir, the discharge radius, thickness, and porosity are identiﬁed as the primary factors inﬂuencing gas well productivity. Furthermore, seven horizontal ﬂow charts illustrating the double-layer gas reservoir and ﬁve horizontal ﬂow charts depicting single-factor variations in the double-layer gas reservoir were constructed. These charts provide a clear visualization of the impact of each reservoir factor on stratiﬁcation’s contribution rate. In contrast to previous studies, this novel approach presents a comprehensive optimization framework that ranks the inﬂuence weights of individual factors and identiﬁes the most signiﬁcant factors impacting multi-layer gas reservoirs. The presented method also serves as a foundation for the subsequent selection of multi-layer gas reservoirs, formulation of gas well stimulation measures, and efﬁcient development.


Introduction
Due to the formation of sedimentary layers in various geological ages, many reservoirs exhibit multiple layer characteristics. The permeability and porosity vary across these layers, which can either be connected or disconnected. In cases where there is strong interlayer leakage flow and good inter-layer connectivity, the dynamic characteristics resemble those observed in a single-layer scenario. However, if the layers are only connected through the wellbore, they clearly demonstrate multi-layer characteristics [1].
The characteristics are significantly different between multi-layer gas reservoirs and single-layer gas reservoirs. This is evident in the disparate depletion rates observed across each production layer, leading to variations in pressure and occurrences of interlayer channeling as well as wellbore backflow [2]. To meet production requirements, many production wells employ a multi-layer combined production approach. During production, the output of multi-layer commingled wells may occasionally be lower than the sum of individual layer outputs, which contradicts expectations and often perplexes reservoir engineers in the field. Investigating seepage characteristics and productivity factors of layered reservoirs is essential for practical analysis and interpretation of multi-layer commingled gas wells, guiding the selection of appropriate production pressure differentials, and facilitating development of multi-layer gas fields.
The investigation of unstable seepage flow in multi-layer commingled mining has been conducted since the 1960s, with numerous valuable studies undertaken by esteemed scholars [3][4][5]. Lefkovits et al. pioneered the application of analytical methods to investigate a multi-layer system model devoid of interlayer fouling, obtaining solutions for bottomhole pressure and production in a circular closed reservoir with an arbitrary number of layers, and conducting a comprehensive analyzing specifically focusing on the two-layer case. While the study holds significant importance, its practical utility is constrained by its failure to consider inter-layer crossflow. Russell and Prats conducted a comprehensive investigation on the multi-layer system model with crossflow, resulting in a novel solution for fixed bottom hole pressure, effectively addressing the limitations of previous research. The study conclusively demonstrated that inter-layer crossflow significantly influences the dynamic response of pressure during the transition stage from no crossflow to a singlelayer homogeneous response. Additionally, scholars have extensively explored diverse initial pressure, boundary, porosity, and permeability conditions by integrating them into various models for further examination [6][7][8]. Papadopulos successfully resolved the issue of disparate pressures and pressure distributions across layers by obtaining an exact solution for a two-layer homogeneous and infinitely bounded aquifer model. Kucuk et al. developed an analytical solution applicable to scenarios involving distinct initial pressures and external boundary conditions in each layer, including complex cases featuring partial perforation, fractured wells, dual media, and external boundary conditions. The aforementioned scholars overlooked the inclusion of the wellbore storage coefficient and skin factor in the study. However, as research progressed, it became evident to many scholars that incorporating these factors into the model is crucial, leading to numerous subsequent investigations [9][10][11][12][13]. Tariq et al. expanded the well-testing model for a multi-layer non-channeling flow system, incorporating both the effects of the wellbore storage effect and skin factor. Larson addressed the challenge of multi-layer commingled production with varying initial conditions, but it was assumed to be under constant rate conditions in the wellbore, and there still exist challenges in dealing with the problem of wellbore storage. Bourdet proposed a two-layer channeling reservoir model that incorporates the effects of the wellbore storage and skin effect. In the subsequent year, Mavor utilized a combinatorial technique to consider the influence of the wellbore storage effect, and applied the concept of parallel resistance to solve the multi-layer connected reservoir model. Rahman and Mattarm introduced a new analytical model that accounts for both the skin factor and wellbore storage effects in order to address instantaneous flow problems associated with commingled production in tight gas reservoirs featuring different initial pressures. This model is suitable for handling variations in production steps before and during production and construction stages, as well as constant pressure or non-flow conditions at the outer boundary of each layer.
Additionally, researchers have studied the calculation model for flow in each layer under varying pressure factors and obtained a series of results during this process [14,15]. Several scholars have utilized the relationship between flow rate at a constant pressure and pressure at a constant flow rate [14,[16][17][18]. By considering the fundamental condition that the flow rates of individual layers are independent from each other during constant pressure production, the respective layer's flow rate is determined, followed by determining the wellbore pressure. Milad et al. established a hydraulic coupling parameter model for multi-layer shale gas reservoirs using an effective iterative numerical simulation method to couple wellbore and gas reservoirs. A methodology is proposed for modeling and simulating commingled production from multiple shale-gas layers while accounting for changes in both pressure and flow within each layer.
Many scholars have conducted extensive research on production performance and prediction models in order to enhance the utilization of multi-layer commingled production technology for high-efficiency development in practical production. Ei-Banbi and Robert developed a stratified steady flow model for multi-layer commingled production data of tight gas reservoirs by combining the material balance equation with the gas steady flow equation [19]. This approach enables analysis and prediction of production performance for multi-layer commingled production, obtaining geological reserves of natural gas (OGIP) for each layer in tight gas reservoirs and evaluating the relative importance of each layer. Building upon Ei-Banbi and Rober's established multi-layer commingled production model, Arevalo-Villagran improved the productivity prediction method for multi-layer commingled gas wells by further incorporating the actual boundary flow equation into the material balance equation, thereby proposing a novel multi-layer commingled production model capable of matching and predicting production outcomes [20]. This model can be used to calculate natural gas geological reserves, productivity index, and production rate of each layer while also enabling predictions regarding the individual layer's future output.
As widely recognized, the analysis of the interlayer interference mechanism in multilayer commingled production systems has been a focal point of research on multi-layer commingled production models for decades. Over the course of several decades, numerous scholars have conducted studies on this issue utilizing numerical simulations and other methods, resulting in a series of valuable findings. After analyzing the main factors contributing to interlayer interference, Zhang et al. concluded that they believed multilayer commingled mining and differences in pressure systems were conditions for such interferences to occur [21]. Wang et al. [22,23] developed a methodology to calculating the interlayer interference coefficient of gas wells based on the binomial productivity formula and wellbore pressure calculation model derived from a physical model of multi-layer commingled production. Over the past decade, numerous scholars have extensively investigated the factors influencing interlayer interference in multi-layer commingled mining through the establishment of diverse mathematical and physical models, as well as utilizing numerical simulations, laboratory tests, and other methodologies. According to variations in research methodologies, the analysis can be categorized into three main approaches: (1) Through the establishment of a two-layer model without crossflow between layers, a physical model of multi-layer commingled production in gas wells, a rate transient analysis (RTA) and a structured RTA model of gas wells, a steady-state two-phase pipe flow model, and a radial multi-layer commingled production numerical model [24][25][26][27][28], the factors influencing interlayer interference were analyzed, and it was determined that the main factors affecting interlayer interference were interlayer pressure difference, reservoir physical properties, and interlayer span. Furthermore, the impact of each layer's properties on multi-layer commingled production was evaluated, including permeability, porosity, thickness, initial pressure, and compressibility. (2) By conducting a visible multi-tube water flooding experiment, Huang et al. successfully quantified the interlayer interference in the multi-layer commingled production of offshore heavy oil reservoirs and established a theoretical model to predict directional well productivity [29]. (3) Zhao et al. employed the grey lattice Boltzmann method to evaluate how permeability and pressure differences affect interlayer interference in multi-layer commingled production [30]. It can be inferred that the interlayer interference coefficient exhibits a linear increase with the rise in pressure differential during multi-layer commingled production of gas reservoirs. Moreover, reservoirs with lower permeability exhibit more pronounced levels of interference.
In recent years, conventional oil production has remained stable, while advancements in the oil industry have significantly boosted unconventional oil production. Compared to conventional oil and gas reservoirs, unconventional oil and gas reservoirs exhibit a more intricate geological structure. Therefore, it is imperative to consider various factors such as reservoir density and heterogeneity when investigating the seepage behavior of these reservoirs. By incorporating the dynamic characteristics of three distinct pore systems, namely inorganic matter, organic matter, and natural fractures, Yan et al. [31] have successfully developed a micro-scale multi-pore model for fluid flow in shale reservoirs. This advancement not only enhances our understanding of the intricate flow mechanism but also paves the way for its potential extension to the reservoir scale through technological advancements. In order to investigate the flow mechanism in the re-reservoir of tight oil reservoirs, Wang and Yan et al. [32,33] conducted a series of studies on modeling the composition of tight oil. They proposed a comprehensive model for tight oil reservoir composition that incorporates the influence of nanopores. This model accurately simulates the rock and fluid characteristics affected by nanopores, resolves inconsistencies in the gas-oil ratio (GOR), significantly simplifies history matching procedures, and enhances recovery prediction reliability.
Numerical simulation methods are crucial for studying oil and gas reservoir seepage behavior. In the past decade, researchers have focused on enhancing these methods through novel algorithms. Tarip, Xu, Gudala et al. [34][35][36][37][38] have made significant advancements using the integrating inversion neural network (INN), particle swarm optimization evolutionary algorithm (PSO), and deep learning (DL) models to investigate reservoir seepage mechanisms and management frameworks. This represents a recent research breakthrough that requires further exploration to provide a more comprehensive description of reservoir characteristics.
In summary, the analysis of interlayer interference mechanism and influencing factors in a multi-layer commingled gas reservoir is typically done through three methods: indoor experiments, numerical simulations, and production model establishment. However, these methods have some drawbacks such as high costs for laboratory experiment simulation and a lack of parameters in production models or quantitative characterization in some approaches. The present paper proposes a novel mathematical model for multi-layer commingled production taking into account multiple influencing factors. By integrating the unstable seepage model of multi-layer commingled production with orthogonal experimental design, a comprehensive analysis is conducted on the proportion of each control factor affecting the productivity of commingled gas reservoirs, leading to a profound understanding of the unique seepage behavior in such reservoirs. In comparison with previous research, this model considers the influence of multiple factors on productivity and employs orthogonal testing to rank their importance in relation to multi-layer gas reservoirs. This provides a foundation of gas well stimulation measures, and efficient development practices for multi-layer gas reservoirs.
The following provides a concise overview of this paper. In Section 2, the governing equations and analytical solutions of the unstable seepage model of multi-layer commingled production are presented. Section 3 conducts an analysis on the combination of orthogonal test design, determining both the primary and secondary order of each influencing factor as well as their respective influence weights. Finally, in Section 4, conclusive remarks are provided. Figure 1 depicts the schematic diagram of a physical model for multi-layer combined gas production. In this model, it is assumed that there exists a combined gas production well at the center of a closed gas reservoir with initial pressure p ij , permeability k j , porosity ϕ j , thickness h j , and drainage radius r ej maintained at constant production rate. Under the general modeling condition, the governing equation for micro-compressible fluids is established:

Establishment of Constant Production Mathematical Model
Water 2023, 15, x FOR PEER REVIEW 5 of 20 Figure 1 depicts the schematic diagram of a physical model for multi-layer combined gas production. In this model, it is assumed that there exists a combined gas production well at the center of a closed gas reservoir with initial pressure pij, permeability kj, porosity φj, thickness hj, and drainage radius rej maintained at constant production rate. Under the general modeling condition, the governing equation for micro-compressible fluids is established:

Establishment of Constant Production Mathematical Model
Initial conditions: Outer boundary conditions: Inner boundary conditions:  Initial conditions: Outer boundary conditions: Inner boundary conditions: where p ji -original formation pressure for layer j, MPa; p pj -gas pseudo pressure for layer j, MPa; ϕ j -porosity for layer j; k j -permeability for layer j, µm 2 ; k fdj -recharge permeability for layer j, µm 2 ; h j -thickness for layer j, m; c tj -comprehensive compression coefficient for layer j, 1/MPa; r ej -drainage radius for layer j, m; r fd -drainage radius for layer j, m; r-wellbore radius, m; p js (t)-bottomhole pressure, MPa; q jsc (t)-surface production for layer j, m 3 /D; and q sc -surface production of a gas well, m 3 /D.

Laplace Solution of Constant Production Equation
The variable flow problem is summarized for the aforementioned Equations (1)- (6). In accordance with the Duhamel superposition principle, the following relations are established for the j th layer: where p pjsr represents the solution to the problem subject to dimensionless flow per unit conditions. Upon applying dimensionless and Laplace transform, we obtain: If the dimensionless flow of a single layer is all 1, that is, q Jd = 1, the solution of the pull-space pressure under fixed production can be obtained through Laplace transform as follows: In the formula: In the case of a closed boundary, λ j = 0, while for a constant pressure boundary, λ j → ∞. I 0 (v) and I 1 (v) are the first class of deformed Bessel functions, whereas K 0 (v) and K 1 (v) belong to the second class of deformed Bessel functions. The calculation of single-layer liquid production q jD (t D ) can be performed using Formula (8).

Real-Time Domain Solution of Constant Production Problem
Define the average pressure as: where V jT -total volume of discharge area, m 3 and V w -discharge shaft volume, m 3 . Then, for the plane radial flow case: Substitute Formula (8) into Formula (13) after Laplace transformation to obtain: Substitute Formula (9) into the above formula and invert by Stehfest to obtain the real-time domain solution of the average pressure: After dimensionalization, it can be obtained: (15) where the left side of the equation represents the average pseudo pressure drop in layer j, while the right side is the cumulative production. Taking a two-layer combined reservoir as an example, the permeability of each layer is maintained at a constant level to investigate the influence of the storage coefficient on layer flow behavior. Figure 2 illustrates the relationship between layer flow and the storage coefficient under different permeability conditions with equal initial pressure. When both layers have equal permeability, the layer with a higher storage coefficient exhibits a greater contribution to flow. Conversely, when the permeability of the two layers is unequal, the layer with a higher permeability and storage coefficient has a higher flow contribution. Additionally, in cases where there is a high permeability but low storage coefficient, there exists an intersection point in the early stage of production on the flow contribution curve for that particular layer; however, during later stages of production, it is observed that layers with a high permeability exhibit higher flow contributions.

Permeability Coefficient
The storage coefficient of each layer remains unchanged to investigate the impact of the permeability coefficient on the layer flow. Figure 3 illustrates the relationship between the layer flow and permeability coefficient for each layer, considering equal initial pressure but different storage coefficients. The influence of the permeability coefficient on layer flow exhibits a similar pattern to that of the storage coefficient. Layers with a high permeability coefficient and storage coefficients contribute more significantly throughout the production process, while layers with a low permeability but high storage coefficients exhibit intersecting contribution curves during production. However, the overall impact of permeability and storage coefficients on layer flow in the entire production process remains unclear.

Dimensionless Initial Pressure
The results depicted in Figure 4 indicate that the initial pressure exhibits heterogeneity, leading to the occurrence of backflow during the early stages of well opening. The magnitude of backflow is directly proportional to both the permeability and pressure difference, while inversely proportional to the storage coefficient. Moreover, a greater disparity in initial pressure intensifies the severity of the backflow phenomenon. However, as mining time progresses, the bottom hole pressure decreases gradually, resulting in a rapid disappearance of backflow.

Effect of Formation Parameters on Layered Contribution Rate
Taking the double-layer gas production well of a tight gas reservoir in the Changqing Oilfield as an example, the average daily gas production rate is 1 × 10 4 m 3 /d. The thickness of the gas reservoir ranges from 2 to 10 m, with permeability varying between 2 and 10 × 10 −3 µm 2 and porosity ranging from 8% to 12%. The initial pressure of the lowpressure gas reservoir is 17.8 MPa, while that of the high-pressure gas reservoir is 18 MPa. Additionally, the drainage radius is 100-500 m. Based on these parameters' distribution within four groups in the gas reservoir, five levels are designed for each group in the low-pressure gas reservoir, as shown in Table 1. greater contribution to flow. Conversely, when the permeability of the two layers is unequal, the layer with a higher permeability and storage coefficient has a higher flow contribution. Additionally, in cases where there is a high permeability but low storage coefficient, there exists an intersection point in the early stage of production on the flow contribution curve for that particular layer; however, during later stages of production, it is observed that layers with a high permeability exhibit higher flow contributions.

Permeability Coefficient
The storage coefficient of each layer remains unchanged to investigate the impact of the permeability coefficient on the layer flow. Figure 3 illustrates the relationship between the layer flow and permeability coefficient for each layer, considering equal initial pressure but different storage coefficients. The influence of the permeability coefficient on   ity, leading to the occurrence of backflow during the early stages of well opening. The magnitude of backflow is directly proportional to both the permeability and pressure difference, while inversely proportional to the storage coefficient. Moreover, a greater disparity in initial pressure intensifies the severity of the backflow phenomenon. However, as mining time progresses, the bottom hole pressure decreases gradually, resulting in a rapid disappearance of backflow.

Effect of Formation Parameters on Layered Contribution Rate
Taking the double-layer gas production well of a tight gas reservoir in the Changqing Oilfield as an example, the average daily gas production rate is 1 × 10 4 m 3 /d. The thickness of the gas reservoir ranges from 2 to 10 m, with permeability varying between 2 and 10 × 10 −3 μm 2 and porosity ranging from 8% to 12%. The initial pressure of the low-pressure gas reservoir is 17.8 MPa, while that of the high-pressure gas reservoir is 18 MPa. Additionally, the drainage radius is 100-500 m. Based on these parameters' distribution within four groups in the gas reservoir, five levels are designed for each group in the low-pressure gas reservoir, as shown in Table 1. In order to investigate the influence of single-factor permeability on production in stratified gas wells, we designed parameters for a low-pressure gas reservoir and established five levels of permeability, as shown in Table 2.  In order to investigate the influence of single-factor permeability on production in stratified gas wells, we designed parameters for a low-pressure gas reservoir and established five levels of permeability, as shown in Table 2.  Figure 5 shows the relationship between the layer flow and permeability at each layer with other parameters held constant. When the permeability of the low-pressure gas layer is lower than that of the high-pressure gas layer, the contribution of layer flow from the high-pressure gas layer surpasses that from the low-pressure gas layer throughout the entire production process. In cases where there exists a high-permeability low-pressure gas layer and a low-permeability high-pressure gas layer, the contribution curve of layer flow from the high-pressure gas layer intersects with that of the early stage in the production process, but the contribution of the layer flow of the high-pressure gas layer is generally higher in the whole production process.

2
Effect of initial pressure on layered contribution rate.
In order to study the influence of initial pressure on production in stratified gas wells, a low-pressure gas reservoir was selected as the experimental setting, with five different levels of initial pressure established as outlined in Table 3.
high-pressure gas layer surpasses that from the low-pressure gas layer throughout the entire production process. In cases where there exists a high-permeability low-pressure gas layer and a low-permeability high-pressure gas layer, the contribution curve of layer flow from the high-pressure gas layer intersects with that of the early stage in the production process, but the contribution of the layer flow of the high-pressure gas layer is generally higher in the whole production process.  Effect of initial pressure on layered contribution rate.
In order to study the influence of initial pressure on production in stratified gas wells, a low-pressure gas reservoir was selected as the experimental setting, with five different levels of initial pressure established as outlined in Table 3. The relationship between the layer flow rate and the initial pressure at each layer, while keeping other parameters constant, is illustrated in Figure 6. With the increase in the initial pressure of the low-pressure gas layer, the pressure difference between the two layers decreases, the contribution of the low-pressure gas layer during early production stages gradually intensifies but diminishes in later stages. In the whole production process, there is generally a significant contribution from the high-pressure gas layer flow. Overall, it can be inferred that variations in the initial pressure of the low-pressure gas layer have a limited impact on its contribution to overall flow.  The relationship between the layer flow rate and the initial pressure at each layer, while keeping other parameters constant, is illustrated in Figure 6. With the increase in the initial pressure of the low-pressure gas layer, the pressure difference between the two layers decreases, the contribution of the low-pressure gas layer during early production stages gradually intensifies but diminishes in later stages. In the whole production process, there is generally a significant contribution from the high-pressure gas layer flow. Overall, it can be inferred that variations in the initial pressure of the low-pressure gas layer have a limited impact on its contribution to overall flow.  3 Effect of porosity on layered contribution rate.
In order to study the influence of single factor porosity on stratified gas well production, the parameters of the low-pressure gas reservoir were designed and five porosity levels were set, as shown in Table 4.  3 Effect of porosity on layered contribution rate.
In order to study the influence of single factor porosity on stratified gas well production, the parameters of the low-pressure gas reservoir were designed and five porosity levels were set, as shown in Table 4.  Figure 7 shows the relationship between the layer flow and porosity for each layer with other parameters held constant. With the increase in the porosity of the low-pressure gas reservoir, the contribution degree of the layer flow increases significantly in the middle production period, and there is little difference between the early and late production period. In general, the contribution degree of bed flow increases with the increase in porosity.

3
Effect of porosity on layered contribution rate.
In order to study the influence of single factor porosity on stratified gas well production, the parameters of the low-pressure gas reservoir were designed and five porosity levels were set, as shown in Table 4.  Figure 7 shows the relationship between the layer flow and porosity for each layer with other parameters held constant. With the increase in the porosity of the low-pressure gas reservoir, the contribution degree of the layer flow increases significantly in the middle production period, and there is little difference between the early and late production period. In general, the contribution degree of bed flow increases with the increase in porosity.

4
Effect of formation thickness on layered contribution rate.

4
Effect of formation thickness on layered contribution rate.
In order to study the influence of single-factor layer thickness on stratified gas well production, the parameters of the low-pressure gas layer were designed, and five horizontal layers were set, as shown in Table 5.  Figure 8 shows the relationship between the layer flow and thickness of each layer with other parameters held constant. When the thickness difference between the two layers is ≥4 m, the high-pressure layer thickness exhibits a higher contribution to layer flow throughout the production process. Additionally, in the early stages of production, there is an intersection between the contribution curves of low-pressure layer thickness and low-pressure gas layer thickness. For a thickness difference of 2-4 m (i.e., when the lowpressure gas layer has a thickness of 4-6 m), the contribution degree of the low-pressure gas layer surpasses that of the high-pressure gas layer. Towards late production periods, both layers tend to have equal contributions with a degree close to 0.5. throughout the production process. Additionally, in the early stages of production, there is an intersection between the contribution curves of low-pressure layer thickness and low-pressure gas layer thickness. For a thickness difference of 2-4 m (i.e., when the lowpressure gas layer has a thickness of 4-6 m), the contribution degree of the low-pressure gas layer surpasses that of the high-pressure gas layer. Towards late production periods, both layers tend to have equal contributions with a degree close to 0.5.

5
Effect of drainage radius on layered contribution rate.
In order to study the influence of single-factor drainage radius on stratified gas well production, the parameters of the low-pressure gas reservoir were designed, and five levels of drainage radius were set, as shown in Table 6.  Figure 9 shows the relationship between the layer flow and drainage radius for each layer with other parameters held constant. As the drainage radius of the low-pressure gas layer increases, the contribution degree of layer flow gradually rises; however, throughout the entire production process, the contribution degree of the high-pressure gas layer flow remains consistently higher than that of the low-pressure gas layer.

5
Effect of drainage radius on layered contribution rate.
In order to study the influence of single-factor drainage radius on stratified gas well production, the parameters of the low-pressure gas reservoir were designed, and five levels of drainage radius were set, as shown in Table 6.  Figure 9 shows the relationship between the layer flow and drainage radius for each layer with other parameters held constant. As the drainage radius of the low-pressure gas layer increases, the contribution degree of layer flow gradually rises; however, throughout the entire production process, the contribution degree of the high-pressure gas layer flow remains consistently higher than that of the low-pressure gas layer.
The orthogonal experimental design is a crucial discipline within statistical mathematics. The fundamental tool of the orthogonal method is the orthogonal table, which possesses two distinctive characteristics: (1) it ensures an equal number of different levels for each participating factor in the experiment and the number of combinations; (2) each

Multi-Factor Sensitivity Analysis 1
Orthogonal experiment design.
The orthogonal experimental design is a crucial discipline within statistical mathematics. The fundamental tool of the orthogonal method is the orthogonal table, which possesses two distinctive characteristics: (1) it ensures an equal number of different levels for each participating factor in the experiment and the number of combinations; (2) each level of any given factor exhibits closely comparable test conditions, enabling an effective comparison among factors. Therefore, employing orthogonal tables for experimental design not only reduces the number of experiments but also enhances their representativeness, thereby meeting the requirements of comprehensive experimentation more effectively.
High-pressure gas reservoir parameters were kept constant, and the contribution ratio of low-pressure gas reservoir production was taken as the test result, as shown in Table 7, to explore the proportion of each influencing factor under the interaction of multiple factors. Subsequently, range analysis and analysis of variance were conducted on these results to determine the primary and secondary order of influencing factors, as well as the influence weight of the single-layer production contribution ratio under multiple influencing factors. The analytical process is illustrated in Figure 10.

2
Range analysis of orthogonal test results.
Through range analysis, a more intuitive comparison can be made between the primary and secondary factors influencing this test. The underlying principle is that variations in the k value, representing the average of i-level data for factor A, are generally attributed to different levels of factor A. The range (R A ) between k values for different levels of factor A can be considered an approximate measure of how changes in factor A levels impact experimental results, with R A positively correlating with the significance of factor A.  2 Range analysis of orthogonal test results.
Through range analysis, a more intuitive comparison can be made between the primary and secondary factors influencing this test. The underlying principle is that variations in the k value, representing the average of i-level data for factor A, are generally attributed to different levels of factor A. The range (RA) between k values for different levels of factor A can be considered an approximate measure of how changes in factor A levels impact experimental results, with RA positively correlating with the significance of factor A.
The range analysis results are presented in Table 8, with the k value representing the average sum of each test index. Among them, the optimal levels are observed for the fifth level of thickness, the fifth level of permeability, the second level of porosity, the first level of initial pressure, and the second level of drainage radius. The range R indicates that the main and second order of influencing factors is drainage radius > permeability > thickness > initial pressure > porosity. Furthermore, the trend chart of contributing factors of gas reservoir production in Figure 11 also reflects the above conclusion. The three most important factors are drainage radius, permeability, and thickness, and porosity and initial pressure have little influence on production. The range analysis results are presented in Table 8, with the k value representing the average sum of each test index. Among them, the optimal levels are observed for the fifth level of thickness, the fifth level of permeability, the second level of porosity, the first level of initial pressure, and the second level of drainage radius. The range R indicates that the main and second order of influencing factors is drainage radius > permeability > thickness > initial pressure > porosity. Furthermore, the trend chart of contributing factors of gas reservoir production in Figure 11 also reflects the above conclusion. The three most important factors are drainage radius, permeability, and thickness, and porosity and initial pressure have little influence on production. The isoline map presented in Figure 12 illustrates the contribution ratio of the two primary factors, namely the drainage radius and permeability, to the gas reservoir production. In this representation, the blue region is the area with a low contribution ratio, and the red area is the area with a high contribution ratio. As depicted in the figure, an increase in both drainage radius and permeability leads to a corresponding elevation in the production contribution ratio for this particular layer. 3 Variance analysis of Orthogonal test results.
Variance analysis is used to infer the significance of the difference in the overall mean value of the corresponding factors at different levels by analyzing the difference between the data fluctuation caused by each factor and the data fluctuation caused by the error. Given the complexity of this experiment involving multiple factors, batch numerical analysis software was employed for efficient processing of variance and range analyses in this study.  Figure 11. Trend chart of factors contributing to gas reservoir production.
The isoline map presented in Figure 12 illustrates the contribution ratio of the two primary factors, namely the drainage radius and permeability, to the gas reservoir production. In this representation, the blue region is the area with a low contribution ratio, and the red area is the area with a high contribution ratio. As depicted in the figure, an increase in both drainage radius and permeability leads to a corresponding elevation in the production contribution ratio for this particular layer.  Variance analysis of Orthogonal test results.
Variance analysis is used to infer the significance of the difference in the overall mean value of the corresponding factors at different levels by analyzing the difference between the data fluctuation caused by each factor and the data fluctuation caused by the error. Given the complexity of this experiment involving multiple factors, batch numerical analysis software was employed for efficient processing of variance and range analyses in this study. Table 9 shows the results of analysis of variance, where the F value represents the level of significance and the p value indicates the magnitude of difference. The findings reveal that the maximum F value of the drainage radius is 3.89, and the significance degree accounts for 38.71%, as shown in Figure 13. It was followed by thickness and permeability, whose F-values were 2.84 and 2.48, respectively, and the significant degree accounted for 28.26% and 24.68%. However, porosity and initial pressure do not exert significant influence on yield contribution. The order of importance for significant effects is as follows:  Table 9 shows the results of analysis of variance, where the F value represents the level of significance and the p value indicates the magnitude of difference. The findings reveal that the maximum F value of the drainage radius is 3.89, and the significance degree accounts for 38.71%, as shown in Figure 13. It was followed by thickness and permeability, whose F-values were 2.84 and 2.48, respectively, and the significant degree accounted for 28.26% and 24.68%. However, porosity and initial pressure do not exert significant influence on yield contribution. The order of importance for significant effects is as follows: drainage radius > thickness > permeability > initial pressure > porosity, and the three most important factors were drainage radius, permeability, and thickness, which was consistent with the results of range analysis.

Conclusions
(1) The paper establishes an unsteady seepage flow model for the combined production of multilayer gas reservoirs in vertical wells. By utilizing the Duhamel convolution principle, the variable flow problem is transformed into a constant flow problem. Subsequently, the Laplace spatial solution is obtained through Laplace transform, and finally, the real time-domain solution is derived using the Stehfest inversion. This comprehensive model takes into account various parameters such as permeability, reservoir thickness, porosity, and initial pressure to fully assess productivity.
(2) By integrating the unstable percolation model of multi-layer gas reservoirs with an orthogonal test design, this study comprehensively considers the proportion of each influencing factor under the combined influence of multiple factors. In the case of the double-layer gas production well of a tight gas reservoir in the Changqing Oilfield in this paper, where there is a small initial pressure difference between the gas reservoirs and no backflow phenomenon was found, it is found that the impact on productivity contribution is minimal due to the small permeability scale and span. The test results were analyzed through a combination of an orthogonal experimental design, range analysis, and variance analysis. The results of the range analysis and variance analysis indicate that the drainage radius has the greatest impact, followed by permeability, thickness, initial pressure, and porosity. The findings demonstrate that, within this gas reservoir, the dominant factors influencing productivity are the drainage radius, permeability, and thickness, whereas initial pressure and porosity exhibit a negligible impact on productivity.
(3) The model proposed in this paper provides valuable practical guidance for current multi-layer combined gas production reservoirs. It effectively orders the main and secondary influencing factors under the common influence of multi-layer gas reservoirs with different characteristics, thus offering theoretical guidance for subsequent joint drilling and abandonment of the secondary productive contribution layer in multi-layer gas reservoirs.

Conclusions
(1) The paper establishes an unsteady seepage flow model for the combined production of multilayer gas reservoirs in vertical wells. By utilizing the Duhamel convolution principle, the variable flow problem is transformed into a constant flow problem. Subsequently, the Laplace spatial solution is obtained through Laplace transform, and finally, the real time-domain solution is derived using the Stehfest inversion. This comprehensive model takes into account various parameters such as permeability, reservoir thickness, porosity, and initial pressure to fully assess productivity.
(2) By integrating the unstable percolation model of multi-layer gas reservoirs with an orthogonal test design, this study comprehensively considers the proportion of each influencing factor under the combined influence of multiple factors. In the case of the double-layer gas production well of a tight gas reservoir in the Changqing Oilfield in this paper, where there is a small initial pressure difference between the gas reservoirs and no backflow phenomenon was found, it is found that the impact on productivity contribution is minimal due to the small permeability scale and span. The test results were analyzed through a combination of an orthogonal experimental design, range analysis, and variance analysis. The results of the range analysis and variance analysis indicate that the drainage radius has the greatest impact, followed by permeability, thickness, initial pressure, and porosity. The findings demonstrate that, within this gas reservoir, the dominant factors influencing productivity are the drainage radius, permeability, and thickness, whereas initial pressure and porosity exhibit a negligible impact on productivity.
(3) The model proposed in this paper provides valuable practical guidance for current multi-layer combined gas production reservoirs. It effectively orders the main and secondary influencing factors under the common influence of multi-layer gas reservoirs with different characteristics, thus offering theoretical guidance for subsequent joint drilling and abandonment of the secondary productive contribution layer in multi-layer gas reservoirs.

Data Availability Statement:
The data presented in this study are available on request from the corresponding authors. The data are not publicly available due to the continuation of a follow-up study by the authors.

Conflicts of Interest:
The authors declare no conflict of interest.