An alternative approach to match field production data from unconventional gas-bearing systems

Nowadays, the unconventional gas-bearing system plays an increasingly important role in energy market. The performances of the current history-matching techniques are not satisfied when applied to such systems. To overcome this shortfall, an alternative approach was developed and applied to investigate production data from an unconventional gas-bearing system. In this approach, the fluid flow curve obtained from the field is the superposition of a series of Gaussian functions. An automatic computing program was developed in the MATLAB, and both gas and water field data collected from a vertical well in the Linxing Block, Ordos Basin, were used to present the data processing technique. In the reservoir study, the automatic computing program was applied to match the production data from a single coal seam, multiple coal seams and multiple vertically stacked reservoirs with favourable fitting results. Compared with previous approaches, the proposed approach yields better results for both gas and water production data and can calculate the contributions from different reservoirs. The start time of the extraction for each gas-containing unit can also be determined. The new approach can be applied to the field data prediction and designation for the well locations and patterns at the reservoir scale.


Introduction
The unconventional gas-bearing system, containing continuous accumulations of unconventional natural gases, has attracted much attention with the depletion of conventional natural gases (Feng et al. 2016;Cui et al. 2019). It is widely found in the Ordos Basin, China, and Cooper Basin, Australia. As two or more gas-bearing reservoirs co-exist in the systems bringing more difficulties, current history-matching approaches do not perform well when applied to such system. This work mainly attempts to overcome these difficulties.
The concept of a gas-bearing system can be dated back to the 1970s by Dow (1974) and was further developed by Masters (1979) and Magoon (1994). In the system, a series of gas reservoirs can be developed continuously in adjacent strata in the same basin. With the exhaustion of conventional reservoirs, more researchers investigated unconventional gas-bearing systems (Collett 2002;Curtis 2002;Law 2002;Schmoker 2002;Shurr and Ridgley 2002). The differences exist between unconventional and conventional gas systems (Ayers 2002). Taking coal seam gas as an example, it is a self-sourcing reservoir and stores hydrocarbons primarily in an adsorbed state and secondarily as free gas (Crosdale et al. 1998;Mastalerz et al. 2004). Consequently, the historymatching and analysis approaches of field data for the unconventional gas-bearing system should be different from that in a conventional system (Su et al. 2005).
History-matching methods of field data can be categorised into a mathematics-based approach and a physicsbased approach (Cui et al. 2018c). The mathematics-based approach establishes the decline curve equation and uses Edited by Yan-Hua Sun curve fitting to match the field data, and it has attracted more attention due to its ease of use. The most popular Arps curves (Arps 1945) were classified into three types based on the decline exponent value: a harmonic decline, an exponential decline and a hyperbolic decline. Considering the difference between the conventional and unconventional gas reservoirs, some improved mathematical methods have been developed in recent years such as the power-law exponential decline method (Ilk et al. 2008a, b), the stretched exponential decline method (Valko 2009;Valko and Lee 2010) and the Duong method (Duong 2011).
On the other hand, proponents of the physics-based approach establish a set of partial differential equations (PDEs) to match the bottom hole pressure curve or the production data curve. There are two categories of the physical approach based on the solution method: an analytical method and a numerical method. For the analytical method, the most popular method in recent years is the tri-linear flow model (Brown et al. 2009;Stalgorova and Mattar 2012), simulating flow within the fracture, flow within the stimulated region and flow within the un-stimulated region. The numerical method is popular with the development of computing speed and can be divided into continuous and discrete models. For the continuous model, the dual-porosity models and multi-porosity models are widely used (Schepers et al. 2009;Cipolla et al. 2010;Fan et al. 2010;Biswas 2011), and in these models, the computational notes on the grids can represent different physical meanings. For the discrete model, the widely used method is discrete fracture network (DFN) where the fractures are directly modelled and the computational notes represent either matrix or fracture (Dershowitz et al. 2010;Wang et al. 2018).
When applied to the unconventional gas-bearing system, these approaches perform poorly: (1) for the mathematicsbased approach, most methods only describe the declining trend (Ilk et al. 2008a, b), while the field data from unconventional gas-bearing systems may exhibit multiple peaks because of the overlapped gas-bearing units (Shen et al. 2018); (2) for the physics-based approaches with an analytical solution, they are used to simplify the reservoir property to search for an analytical solution while ignoring the couplings between the gas flow and solid deformation (Ozkan et al. 2011); (3) for the physics-based approaches with a numerical solution, they are usually time-and resource-consuming (Biswas 2011) and the required data and information are not available for all reservoirs (Cipolla et al. 2010), especially for unconventional gas-bearing systems as they contain two or more gas-bearing systems. Furthermore, within the unconventional gas-bearing strata, gas-bearing reservoirs are stacked vertically with a lack of the fluid connectivity between different stratigraphic features (Wang et al. 2015;Shen et al. 2016), causing difficulty in determining the start time of the extraction for each reservoir. Moreover, each gas-bearing unit has its own reservoir pressure, gas transport properties and hydrodynamic conditions (Ebanks 1987;Cui et al. 2018b). Therefore, it is difficult to find an applicable method for all the reservoirs in the same basin. All of factors aforementioned make the matching of the unconventional gas-bearing system difficult.
As stated previously, the current history-matching approaches are insufficient for unconventional gas-bearing systems. To cope with the defects of the current approaches, a mathematical model was proposed with an auto-computing program to match the field data from the unconventional gas-bearing system. The field data collected from the Linxing Block, Ordos Basin, were applied to illustrate the data processing process. Then the automatic program was applied to the reservoir study to verify its applicability. The paper is organised as follows: the mathematical model and autocomputing program are introduced in Sect. 2, the verified case is introduced in Sect. 3, and the applications of the proposed model are introduced in Sect. 4, followed by the conclusions. Figure 1 illustrates the different types of unconventional resources including gas reservoir (coalbed methane, shale gas and tight reservoirs), oil reservoir (oilshale, shale oil and heavy oil) and other reservoirs (gas hydrates). In some sedimentary basins, it may contain one or more specific reservoirs which are stacked vertically. For an unconventional gas-bearing system, it contains at least two gas-bearing units. For the whole strata, each gas-bearing system has unique gas storage and flow properties and can be treated as a separate domain.

Mathematical model
The surface vertical well is usually applied to the coexploitation of an unconventional gas-bearing system as it can reach more gas-containing reservoirs (Li et al. 2018). Usually, a water-based fracturing liquid is used for hydraulic fracturing (Cui et al. 2018b) and gas would be released from the top reservoir to the lower reservoir due to the water flowback process (Cui et al. 2018a).
To present the mathematical model in this work, we established a simple sandwich stratum with two gas-containing reservoirs (tight sandstone gas and coalbed methane) separated by a low-permeability mudstone as shown in Fig. 2. In the figure, the sandstone, mudstone stone and coal seam are vertically stacked with a vertical well going through in the middle. As introduced above, the gases stored in different gas reservoirs release at different times. The light green parts in the well represent the gas mass from different reservoirs. At time t 0 , the gas masses of sandstone and coalbed seam are s 0 and d 0 . At time t, the gas mass of sandstone is represented by s t , while the gas mass of coalbed seam is separated into d tf and d tm because of the different gas transport abilities of fracture and matrix systems in which the subscripts f and m represent the matrix and fracture system. Taking the gas depletion process in sandstone as an example, at the initial state the gas in a sandstone (shown as s 0 in the figure) having a length L s flows in the vertical well with the radius of R and the flow length is L r at the time t.
Taking the gas depletion process in sandstone as an example, the general equation of flow process of the gas mass in the long tube was firstly investigated by Taylor (1954) and can be written as: where the D m is the gas conductivity coefficient and has the same units as the diffusion coefficient, u denotes the flow velocity, and x and r represent the coordinate axis (mass flow direction and well axial direction). For the case where the gas mass forms a very small slug of length L s and the flow path L r is relatively long, the equation can be written as (Kolev 1995): where denotes the dimensionless time as ut/V t , V t represents the characteristic mass per unit volume, X denotes the dimensionless flow distance as x/L, L represents the characteristic length, S represents the dimensionless concentration as sV t /V 0 s 0 and V 0 represents the initial volume of s 0 . The details of their definitions can be found in Kolev's work (1995), and Pe = uLr∕D m is the Peclet number which characterises the dispersion properties of the flow system. The general solution (Kolev 1990) is: where l is the dimensionless length of the initial gas mass. If the gas mass is very small compared to the well volume, it can be approximated by an ideal delta function input. Under the such condition, the solution can be written as:  Fig. 2 Gas depletion process of the unconventional gas-bearing system. R is the inner radius of the well, L s is the length of the mass of sandstone tight gas, and L r represents the length of the flow path. At the initial state (t = 0), the gas masses are in their initial states (s 0 and d 0 ) In fact, the dimensionless concentration (S) represents the mass (gas or water) flow rate measured at the well and Eq. (4) can be written in the general form: where q(t) represents the mass (gas or water) flow rate variations with time. Equation (5) is the standard form of Gaussian function whose graph is a characteristic symmetric bellshaped curve. a, b and c are the curve parameters controlling the shape of the Gaussian function. a is the height of the curve's peak, b is the time of the centre of the peak, and c controls the width of the curve.
As discussed above, the gas flow in the well can be simplified as the dispersion problem in the long tube and its solution is Gaussian function under the assumption that gas mass is very small compared to the well volume. Additionally, the Gaussian functions are widely applied in statistics to describe the normal distributions, in signal processing to define Gaussian filters, in image processing where two-dimensional Gaussians are used for a Gaussian blur and in mathematics to solve heat equations and diffusion equations. In the unconventional gas-bearing system, several gas-bearing units coexist. Furthermore, in some reservoirs, several gas-containing porosity systems coexist and exhibit huge differences in gas transport abilities (Cui et al. 2020). These two factors would lead to the coexistence of the Gaussian function (Cui et al. 2018c). In the full-scale production process, the gas flow rate infield data are considered to be a superposition of several components which can be described by: in which q(t) represents the gas/water flow rate of the whole unconventional gas-bearing system and n denotes the number of Gaussian functions.

Development of the mathematical tool
In the previous section, we have proved that the field data from an unconventional gas-bearing system can be represented by the superposition of a series of Gaussian functions and the Gaussian functions may overlap with each other. In this section, a mathematical tool was proposed to distinguish the overlapped Gaussian function and calculate the exact expressions. The development of the mathematical tool is illustrated in Fig. 3, and the details are introduced below. Savitzky and Golay (1964) proposed a method for data denoising with a polynomial structure. The general equation of the Savitzky-Golay method can be given as follows:

Noise cancellation
where q is the original field data value, q* refers to the resultant field data value, C r is the coefficient for the field data value of the smoothing window and N is the number of convoluting integers and is equal to the smoothing window size 2m + 1. The index j is the running index of the original ordinate data table. The smoothing array consists of 2m + 1 points, where m is the half width of the smoothing window. The coefficients of the Savitzky-Golay equation (C r ) can be obtained directly from Steinier et al. (1972) or calculated based on the equations presented by Madden (1978) and are given in "Appendix 1".

Regression algorithm
Due to the co-existence of the gas-containing reservoir, the field data from the unconventional gas-bearing system usually largely fluctuate, making the data processing difficult. In that situation, data regression is helpful. In this work, polynomial fitting is applied as shown in Eq. (8): where i represents the exponent of the polynomial, p denotes the coefficient of the polynomial and k is the index number.
The least-squares method is applied to find the value of p i (i = 0, 1, 2, …, N) as given in Eq. (9):  The partial differential of p i can be expressed as: Equation (10) can be expanded to give: or expressed as: This gives an expression for p i (i = 0, 1, 2, …, N):

Derivative of the field data
The data derivative is helpful for the case where the important information is hard to extract from the original data. In this work, the data derivative value (especially the second derivative) is the basis of peak detection. The direct regularisation technique was proposed by Roy (2015), and it selects a posteriori optimal regularisation parameter to solve the inverse problem for estimation of first-order derivative spectra. The higher-order derivative spectra can be obtained using the algorithm in the sequel.
The noise contaminated dataset is denoted as q(t) δ , with noise δ within an interval [a, b], and the first derivative is denoted as w(t) which is a continuous function.
The w(t) can be obtained as follows: where ̂ is the first-order differential on w; superscript T represents the transposition of the matrix; A represents the grid interval matrix which can be termed as design matrix (Kreyszig 2011); for the direct method, Aw = q δ ·η > 1 is suitably chosen scalar quantity satisfying an obvious constraint; the parameter η is widely known as the regularisation parameter. The algorithm is as follows.
Note that on exiting the iterative loop while satisfying in the MATLAB environment is given in "Appendix 2".

Peak detection
As mentioned above, the field data from an unconventional gas-bearing system are a combination of a series of Gaussian functions. The method of accessing the rough initial values for iteration is introduced in this section, and the final value is determined by way of a curve fitting algorithm. The initial iteration values can be obtained from the original field data and its second derivative. The details of the approach are illustrated below, and the overlap of two Gaussian functions is taken as an example.
(14) T +̂ T̂ w − T q = 0 Table 1 lists the calculations needed to obtain the initial values and the upper/lower bounds for each parameter. Figure 4 demonstrates the values of b 0 , b l , b u , a 1 and a 2 . a 1 is the value of the input data at x = b 0 . b 0 is calculated by calculating the minimum of the second derivative, b l and b u are the times where the second derivative is zero, and a 2 is the value of the second derivative at x = b 0 .

Curve fitting
To identify the intrinsic parameters a i , b i and c i from the production data, the Levenberg-Marquardt method is introduced. The Levenberg-Marquardt method was proposed by Levenberg (1944) and Marquardt (1963), and it takes advantage of two methods based on different orders of the gradient: "steepest descent" and "Gauss-Newton" (Lampton 1997;Lourakis 2005).
Considering the system of nonlinear equations in this work, we obtain: where the F(t) is the objective function. The iteration for the LM algorithm gives where k is a non-negative regularised parameter, I is the unit matrix and J(t k ) is the Jacobian of F(t) at t k .
3 Verification of the mathematical tool

Geological conditions of the Ordos Basin and the Linxing Block
Ordos Basin is located in the Northwest China and surrounded by the Yellow River on its eastern, northern and western sides (Li et al. 2008). More than fifty types of mineral resources are found in the basin among which the coal, gas and sandstone-hosted uranium are particularly important for national strategic development. The proven coal and gas reserves are 213.3 billion tons and 1789.3 billion cubic metres, respectively, which account for 13.7% and 17.3% of the national reserve (Feng et al. 2016).  The Linxing region is an economically significant area located in the eastern Ordos Basin near the city of Yulin (Fig. 5). A significant volume of unconventional gas (coalbed methane, shale gas and tight sandstone gas) has been found within this region (Ju et al. 2017). After a long and complex multi-cycle geological history, a multi-period and multilayer sandstone reservoir was formed (Fig. 6) (Li et al. 2016). The Taiyuan and Shanxi Formations are located directly up the Benxi Formation. The Lower Shihezi Formation was deposited directly onto the Shanxi Formation along with fluvial sandstones. The Benxi Formation is one of the main coal-bearing strata containing a one-to four-layer coal seam, and the 8# and 9# coal seams are the main target. The Benxi Formation can be divided into Ben 1 Member and Ben 2 Member. The Taiyuan Formation consists of dark grey-black carbonaceous mudstone, black oil light grey medium-fine sandstone and siltstone. The Taiyuan Formation can be divided into Tai 1 Member and Tai 2 Member. The Shanxi Formation consists of two to five coal layers, and the 4# and 5# coal seams are the main target. The Shanxi Formation can be divided into Shan 1 Member and Shan 2 Member. The Lower Shihezi Formation contains grey mudstone and light grey medium-fine sandstone and can be divided into four members: He 5 to He 8. The Upper Shihezi Formation also contains grey-green mudstone and light fine sandstone and can be separated into He 1 to He 4 Members. Shiqianfeng Formation, 173-286 m thick, is composed of purplish red pebbly coarse sandstone embedded with purplish red to lime green sandy mudstone. The formation can be separated into Qian 1 to Qian 5 Members. The separations of these formations are given in Table 2 from top to bottom (Li et al. 2016).

A verified case study
The production data collected from a vertical well, located in the north part of the Linxing Block, were used for model verification (Li et al. 2018). The well is about 2100 m deep and goes through five gas-bearing strata: three tight gas sandstone layers (Qian 5, He 2 and Tai 2 from the top) and two coal seams (8# and 9# coal seams in Benxi Formation). The gas and water production data with 350 days are displayed in Fig. 7. As shown in the figure, the early gas and water production data oscillate significantly because of the well shutdown; therefore, only long-term production data were used for verification. It is noted that there is only one porosity system in a tight gas sandstone reservoir, while two porosity systems are present in the coalbed methane reservoir (matrix system and a fracture systems); therefore, the total number of Gaussian functions is seven in this case: three for the three tight gas sandstone layers and four for the two coal seams.

Noise cancellation and data regression
Before curve fitting, the Savitzky-Golay method and polynomial fitting were used for noise cancellation and data regression, respectively. The results after using the Savitzky-Golay method and polynomial fitting are plotted in Fig. 8. As illustrated in the figure, the noise cancellation and data regression can reduce data oscillation such as the  rapid decrease between 150 and 160 days and the well shutdown between 240 and 250 days.

Curve fitting
After noise cancellation and data regression, the initial values for iteration were obtained and the Levenberg-Marquardt method was adopted for curve fitting. The results are shown in Fig. 9. Good results were obtained for both gas flow rate and water flow rate (goodness of fit 86% and 80%, respectively).

Discussion
For further discussion, the gas flow rates from different gas reservoirs and porosity systems are illustrated (Fig. 10a). The proportion variations with time are also plotted (Fig. 10b). As shown in Fig. 10a, the gases in different storage systems flow out in sequence: (1) the tight gas in the sandstone flows out first; then, the gas in the coal seam flows out later; (2) the gases stored in the fracture system in the coal seam are first depleted and the gases absorbed in the matrix system diffuse out later. The discrepancy between different gas reservoirs is caused by the water flowback process as the reservoir in the upper strata reaches the gas extraction pressure earlier than the lower reservoir. The difference between different porosity systems arises because of the discrepancy of gas transport abilities as the fracture permeability is usually thousands of times than that of the matrix system. Similar observations can also be made from Fig. 10b: the gas in the sandstone occupies a bigger proportion early in the period, while the coalbed methane does so in the long term. For the time period, we noticed that the gas flow in coal fractures has a shorter period than that in the matrix system because of its larger permeability. To our surprise, the gas in Tai 2 Member has a relatively large period. That is maybe because Tai 2 Member is served as the overlying strata for 8# and 9# coal seams, and most gases in coal seams are flowing into and accumulated in Tai 2 Member. Similarly, the water flow rate and proportion variation with time from different reservoirs and porosity systems are illustrated in Fig. 11. Also, similar conclusions can be obtained: the water in sandstone flows out first for the different gas reservoirs; and the water stored in the fractures is first depleted from the coal seam. The water flow in coal fractures has a shorter timescale than that in the matrix system. It should be noticed that the water in He 2 Member has the longest period from start to the end. The reason for this is that it may include water sourced from other strata. The effect of the number of water sources will be discussed below.

Applications
In this section, the verified mathematical model was applied in field practice. In particular, three cases are introduced: both the gas and water production rate from a single coal seam, both the gas and water production rate from multiple coal seams, and the gas production rate from co-production practice of tight gas and coalbed methane.

Production data from a single coal seam
In this section, the proposed mathematical tool was applied to match the CBM production data from a vertical well in the south of Shanxi Province, China (Ma et al. 2017). Initially, approximately 320 m 3 of pressurised water was injected into the coal seam for hydraulic fracturing to improve the    connectivity to the well and increase the permeability of the reservoir. Both the methane and water were produced through the vertical well. Field history data of methane and water production, from 4 October 2007 to 3 September 2008, were used for model calibration herein. The gas production rate was matched, and the results are shown in Fig. 12. For comparisons, the results of the physics-based approach (Ma et al. 2017), the mathematics-based approach, traditional power-law exponential decline method (Ilk et al. 2008a) and the Duong method (Duong 2011) are drawn. For the mathematics-based approaches, they can well match the early production data, but do not perform well for the later peak. The physics-based approach can match the trend-first decreases and then increases, but its results are unsatisfied. Compared with the previous approaches, our approach is able to well match the results. As shown in the figure, the gas from coal fractures plays the dominant role in the early period, while the gas from the coal matrix takes the full proportion overall in the long-term production process. Similarly, water production data were analysed; as discussed above, the number of water sources plays an important role in the matching results. Firstly, we assume that the water comes only from the coal matrix and coal fractures (two Gaussian functions), then the hydraulic fracturing water is considered (three Gaussian functions), and finally, the water coming from other strata is added (four Gaussian functions). The fitting results are shown in Fig. 13a, and the goodness of fitting is 0.51, 0.67 and 0.75, respectively. Also, the results obtained from Ma's work are also drawn, having a goodness of fitting of 0.45. As shown in the figure, the fitting results are not good when using two or three Gaussian functions, especially for long-term production data, while the use of four Gaussian functions can fit the field data well. The proportions of different resources with four Gaussian     Fig. 12 History-matching results of gas production data from one coal seam: a history-matching results and b proportion with time functions are also plotted in Fig. 13b. As illustrated in the figure, the water from fracturing liquid and coal matrix plays less important roles, while most water comes from the coal fractures and other water-containing strata resources in this case.

Production data from multiple coal seams
In this section, the proposed mathematical tool was applied to match the production data from multiple coal seams. Both gas and water production data were collected from a vertical well in the Ordos Basin, China (Nie et al. 2018). The 5# coal seam in the Shanxi Formation and 8# coal seam in the Taiyuan Formation are its target reservoir, and the 5# coal seam is located above the 8# coal seam. Field history data of methane and water production from 2011 to 2016 were used for model calibration in this study. Firstly, the matching of gas production rate was made and the results are shown in Fig. 14. Two remarkable characteristics can be observed: (1) the gas stored in the 5# coal seam is first depleted, while gas in the 8# coal seam is produced later because of the water flowback process; (2) similar to the above discussion, the gas in the fracture system flows out earlier than that in the matrix system. The water production data were also analysed. The effects of numbers of Gaussian functions were also investigated, and the results are illustrated in Fig. 15a. Four Gaussian functions represent the matrix and fracture systems in both coal seams, the hydraulic fracturing water is considered in five Gaussian functions and the resources from other strata are considered in six Gaussian functions. As shown in the figure, both the set of five and that of six Gaussian functions can fit the field data well, but a poorer result was obtained when using the set of four Gaussian functions. The proportions from different sources are also presented in Fig. 15b. Three characteristics can be inferred from the figure: (1) the hydraulic fracturing water depletes first as characterised by the high water rate; (2) for the coal seam, most water is

Production data from co-production of coalbed methane and tight gas
In this section, the co-production data of CBM and tight gas (TG) from the Linxing Block were applied to the historymatching process (Shen et al. 2018

3
The co-production gas data collected from the six vertical wells in the Linxing Block were applied to history matching. Wells A and B were drilled for the co-production of tight gas in the Shihezi Formation, and both wells are targeting two gas-containing reservoirs. Therefore, the number of Gaussian functions is two for both cases. The history-matching results and the contributions from different gas reservoirs are demonstrated in Fig. 16a, b. As shown in the figure, good fitting results were obtained for both cases.
Wells C and D were drilled for the co-exploitation of coalbed methane and tight gas, and both wells are available to at least three gas-containing gas reservoirs. The fitting results are shown in Fig. 16c, d, and the contributions from different gas reservoirs were also drawn. A good fitting result was obtained in Fig. 16c, while a relatively poor fitting result is shown in Fig. 16d. This is mainly attributed to oscillations in the field data from Well D. On the other hand, the Shan 2 Formation is close to the Shan 1 Formation and the gas that was stored therein may have been depleted.
In the first four cases, the number of Gaussian functions is assumed to be the same as the number of gas-containing systems and good fitting results were acquired. This approach may not always work well. Wells E and G were drilled for the production of tight gas in the Tai 2 Member, and both production data exhibit fluctuations. When one Gaussian function was applied, poor fitting results were obtained for both cases as shown in Fig. 17a, c. As introduced above, the Taiyuan Formation consists of dark grey-black carbonaceous mudstone, black oil light grey medium-fine sandstone and siltstone. In such conditions, we assume that the gas stored in the upper formation (Tai 1 and Shan 2 Members) and lower formation (Ben 2 Member) might be depleted due to the high permeability of the Tai 2 Member. Therefore, the number of Gaussian functions is four (as follows) for the gas in the Tai 2 Member, gas in the Tai 1 and Shan 2 Members, gas in the Ben 2 fracture system and the matrix system. Good history-matching results were obtained with four Gaussian functions as shown in Fig. 17a, c. The contributions from different gas reservoirs for both cases are illustrated in Fig. 17b, d. The gas stored in the Tai 2 Member accounts for a significant proportion in the early period, while gas from the adjacent formation takes the full proportion later. The proportion in the upper formation (Tai 1 and Shan 2 Members) is the biggest. For Well G, the He 7 and Tai 2 Members are its target formation, while when two Gaussian functions were applied poor fitting results were acquired. Then we assume that there would be four Gaussian functions: two for the gas-containing reservoir itself (He 7 Member and Tai 2 Member) and two for its adjacent reservoirs. A good fitting result was obtained as shown in Fig. 17e. Also the contributions from different gas-containing reservoirs are illustrated and most gas comes from the adjacent reservoir as illustrated in Fig. 17f.

Conclusions
In the unconventional gas-bearing system, the varying gasbearing reservoirs are stacked vertically with a lack of fluid connectivity between different stratigraphic features and each gas-bearing unit has its own reservoir properties. Furthermore, traditional history-matching approaches including a mathematics-based approach and physics-based approach do not play well when applied to an unconventional gasbearing system. In this work, an alternative approach was developed and successfully applied at reservoir scale with an automatic calculation program. Based on the calculated results, the main conclusions can be summarised as follows: The production data in a specific gas-containing unit can be described by a Gaussian function since the mass volume is relative small compared with the length of mass flow. Therefore, the coproduction data from the unconventional gas-bearing system were controlled by the supervision of a series of Gaussian functions whose number can be determined based on the amount of the gas-containing units. The adjacent reservoir should be considered when determining the number of Gaussian functions to be used.
Compared with the previously proposed approaches, the proposed method yields better results for both gas and water production data and it can calculate the contributions and proportions from different reservoirs, and determine the star time of the extraction for each gas-containing unit.
The production data characteristic of the unconventional gas-bearing system was investigated using the proposed approach. The flow rate in the early stage is driven mainly from the upper reservoir, while the later flow rate is mainly driven from the lower reservoir because of the water flowback process. Most depleted gases come from the reservoirs where the well is located. When these gases are depleted, the gas stored in the adjacent reservoir would be desorbed and flow into the well.
This work provides an alternative approach to match the gas and water production data from single reservoir and the unconventional gas-bearing system. This approach shows better results and is easy to use; thus, it can be applied to the field data prediction and designation for the well locations and patterns at the reservoir scale. as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.