## Appendix A Derivation of the Mathematical Model  
Based on the mass conservation principle, Darcy's law and real gas equation of state, the governing equations under isothermal circumstance can be obtained respectively. In the inner region, flow conductivity of fracture is greatly higher than matrix, and the double porosity/single permeability model is employed. Natural fracture:  
\[\frac{1}{r}\frac{\partial}{\partial r}\bigg{(}r\frac{\partial p_{pf}}{\partial r }\bigg{)}+\frac{\partial}{\partial z}\bigg{(}\frac{k_{fv}}{k_{fh}}\frac{ \partial p_{pf}}{\partial z}\bigg{)}=\frac{q\mu_{g}c_{t}}{\alpha_{t}k_{fh}} \frac{\partial p_{pf}}{\partial t}-\alpha\frac{K_{m1}}{K_{f}}\frac{K_{f}}{k_{fh }}(p_{pm}-p_{pf}) \tag{1}\]  
Matrix:  
\[\alpha\frac{K_{m1}}{K_{f}}(p_{pm}-p_{pf})+\frac{(qc_{t})_{m}}{(qc_{t})_{f}} \frac{k_{fh}}{K_{f}}\frac{\mu_{g}(\varphi c_{t})_{f}}{k_{fh}}\frac{\partial p_ {pm}}{\partial t}=0 \tag{2}\]  
where, \(K_{m1}\) and \(K_{f}\) are equivalent permeability for matrix and fracture in inner region. As a consequence of anisotropy, the equivalent permeability defined in the following equation can reflect the flowing capacity for fluid mass transfer between matrix and fracture comprehensively. \[K_{m1}=\sqrt{k_{m1v}k_{m1h}},\ K_{f}=\sqrt{k_{fp}k_{fh}} \tag{3}\]  
The pseudopressure functions are defined as:  
\[p_{pf}=\frac{\mu_{g}Z_{i}}{k_{fh}p_{i}}\int\limits_{p_{ac}}^{p_{f}}\frac{k_{ fh}p}{\mu_{g}Z}dp,\ p_{pm1}=\frac{\mu_{g}Z_{i}}{k_{fh}p_{i}}\int\limits_{p_{ ac}}^{p_{m1}}\frac{k_{fh}p}{\mu_{g}Z}dp. \tag{4}\]  
In Equations (1) and (2), the viscosity, compressibility of gas, and the permeability of fracture change with reservoir pressure, hence change with time. Thus, the control equations are still nonlinear. To linearize the control equations, the pseudotime factor can be defined as:  
\[\beta=\frac{1}{t}\frac{\phi_{i}c_{t}\mu_{g}}{k_{fhi}}\int\limits_{0}^{t}\frac {k_{fh}}{\phi c_{t}\mu_{g}}dt \tag{5}\]  
To simplify the governing equations during the derivation, we assume that the following approximate equation is tenable. \[\frac{\left(\phi c_{t}\right)_{m}}{\left(\phi c_{t}\right)_{f}}\approx\frac{ \left(\phi c_{t}\right)_{mi}}{\left(\phi c_{t}\right)_{fi}} \tag{6}\]Based on Equation (A5), the control equations can be further transformed as following:  
\[\alpha\frac{K_{m1}}{K_{f}}(p_{pm1}-p_{pf})+\sqrt{\kappa_{1}}\frac{(\phi c_{t})_{ m1i}}{(\phi c_{t})_{fi}}\frac{\mu_{gi}(\phi c_{t})_{fi}}{k_{fhi}}\frac{\partial p_{pm1}} {\partial(\beta t)}=0\] (A7)  
\[\alpha\frac{K_{m1}}{K_{f}}(p_{pm1}-p_{pf})+\sqrt{\kappa_{1}}\frac{\mu_{gi}( \phi c_{t})_{m1i}}{k_{fhi}}\frac{\partial p_{pm1}}{\partial(\beta t)}=0\] (A8)  
where, \(\kappa_{1}\) is the permeability anisotropy factor of fracture, which is defined in Table A1. In the outer region, the single porosity model is employed. Thus, the control equation for outer region is as following:  
\[\frac{1}{r}\frac{\partial}{\partial r}\bigg{(}r\frac{k_{m2h}}{k_{fhi}}\frac{k_ {fhi}p}{\mu_{g}Z}\frac{\partial p_{m2}}{\partial r}\bigg{)}+\frac{\partial}{ \partial z}\bigg{(}\frac{k_{m2v}}{k_{fhi}}\frac{k_{fhi}p}{\mu_{g}Z}\frac{ \partial p_{m2}}{\partial z}\bigg{)}=\frac{\phi\mu_{g}c_{t}}{\alpha_{t}k_{fhi }}\frac{k_{fhi}p}{\mu_{g}Z}\frac{\partial p_{m2}}{\partial t}\] (A9)  
where,  
\[p_{pm2}=\frac{\mu_{gi}z_{i}}{k_{fhi}p_{i}}\sqrt[p_{m2}]{k_{fhi}p_{i}}dp\] (A10)  
The analogous assumptions are made as:  
\[\alpha\frac{K_{m1}}{K_{f}}(p_{pm}-p_{pf})+\frac{(qc_{t})_{m}}{(qc_{t})_{f}} \frac{k_{fh}}{K_{f}}\frac{\mu_{g}(\varphi c_{t})_{f}}{k_{fh}}\frac{\partial p_ {pm}}{\partial t}=0 \tag{2}\]0 (A11)  
Thus, the control equation for outer region can be written as:  
\[\alpha\frac{K_{m1}}{K_{f}}(p_{pm}-p_{pf})+\frac{(qc_{t})_{m}}{(qc_{t})_{f}} \frac{k_{fh}}{K_{f}}\frac{\mu_{g}(\varphi c_{t})_{f}}{k_{fh}}\frac{\partial p_ {pm}}{\partial t}=0 \tag{2}\]1 (A12)  
By using the dimensionless variables listed in Table A1, the governing equations can be transformed into dimensionless form as follows.