## 3 Mathematical Model
### Coupled Discrete Model  
The flow model for gas reservoirs can associate with the flow model for hydraulic fractures by the following relationship:  
\[\psi_{\rm fD}(r_{\rm D},\theta_{\rm ff},t_{\rm D})=\psi_{\rm D}(r_{\rm D}, \theta_{\rm ff},t_{\rm D}),\ (0<r_{\rm D}<L_{\rm fD}). \tag{25}\]  
Substituting Equations (11) and (A37) into Equation (25), one can derive  
\[\xi_{\rm fD}(r_{\rm D},\theta_{\rm ff},t_{\rm D})=\xi_{\rm D}(r_{\rm D}, \theta_{\rm ff},t_{\rm D}). \tag{26}\]  
The zero-order perturbation solution of the \(\xi_{\rm D}\) is accurate enough in practice, and thus we can obtain  
\[\xi_{\rm D}(r_{\rm D},\theta_{\rm ff},t_{\rm D})=\xi_{\rm D0}(r_{\rm D}, \theta_{\rm ff},t_{\rm D}). \tag{27}\]  
Combining Equations (26) and (27) yields that  
\[\xi_{\rm fD}(r_{\rm D},\theta_{\rm ff},t_{\rm D})=\xi_{\rm D0}(r_{\rm D}, \theta_{\rm ff},t_{\rm D}). \tag{28}\]  
Substituting Equation (28) into Equation (23) leads to  
\[\xi_{\rm wHD}(t_{\rm D})-\xi_{\rm D0}(r_{\rm D},\theta_{\rm ff},t_{\rm D})= \frac{2\pi}{C_{\rm fD}}\bigg{[}r_{\rm D}\int_{0}^{L_{\rm fD}}q_{fD}(r_{\rm D}, \theta_{\rm ff},t_{\rm D}){\rm d}r_{\rm D}-\int_{0}^{r_{\rm D}}\int_{0}^{x^{ \prime}}q_{fD}(x^{\prime\prime},\theta_{\rm ff},t_{\rm D}){\rm d}x^{\prime }{\rm d}x^{\prime}\bigg{]}. \tag{29}\]  
Taking the Laplace transform of Equation (29), one can obtain  
\[\overline{\xi}_{\rm wHD}(s)-\overline{\xi}_{\rm D0}(r_{\rm D},\theta_{\rm ff},s)=\frac{2\pi}{C_{\rm fD}}\bigg{[}r_{\rm D}\int_{0}^{L_{\rm fD}}\overline{q} _{fD}(r_{\rm D},\theta_{\rm ff},s){\rm d}r_{\rm D}-\int_{0}^{r_{\rm D}}\int_{0 }^{x^{\prime}}\overline{q}_{fD}(x^{\prime\prime},\theta_{\rm ff},s){\rm d}x^{ \prime\prime}{\rm d}x^{\prime}\bigg{]}. \tag{30}\]  
Substituting Equation (12) into Equation (30), one can derive  
\[ \begin{array}{c}\overline{\xi}_{\rm wHD}(s)-\sum\limits_{j=1}^{n}\Big{[}\int _{0}^{L_{\rm fD}}\overline{q}_{fD}(\kappa,\theta_{\rm ff},s){\rm K}_{0}(\sqrt{s }\sqrt{r_{\rm D}^{2}+\alpha^{2}-2r_{\rm D}\alpha\cos\left(\theta_{\rm ff}- \theta_{\rm ff}\right)}\big{)}{\rm d}\alpha\Big{]}\\ \qquad\qquad\qquad=\frac{2\pi}{C_{\rm fD}}\bigg{[}r_{\rm D}\int_{0}^{L_{\rm fD }}\overline{q}_{fD}(r_{\rm D},\theta_{\rm ff},s){\rm d}r_{\rm D}-\int_{0}^{r_{ \rm D}}\int_{0}^{x^{\prime}}\overline{q}_{fD}(x^{\prime\prime},\theta_{\rm ff },s){\rm d}x^{\prime\prime}{\rm d}x^{\prime}\bigg{]}.\end{array}  \tag{31}\]  
Employing the Laplace transform of Equations (19) and (20), one can obtain  
\[\sum\limits_{i=1}^{n}\int_{0}^{L_{\rm fD}}\overline{q}_{fD}(r_{\rm D},\theta_{ \rm ff},s){\rm d}r_{\rm D}=\frac{1}{s}. \tag{32}\]  
The coupled model consists of Equations (31) and (32), which can be solved by the numerical discrete method [40]. Each hydraulic fracture is discretized into some segments, and the flow rate in each segment is considered to stay the same at a certain time. The discrete schematic of a hydraulic fracture of a vertical fractured well is shown in Figure 4.  
And then, the discrete forms of Equations (31) and (32) can be obtained as follows  
\[ \begin{split}&\overline{\xi}_{\text{wHD}}(s)-\sum\limits_{i=1}^{n} \sum\limits_{j=1}^{N_{i}}\left[\overline{q}_{\text{fi,jD}}(s)\int_{r_{i,j} \text{D}}^{r_{i,j+10}}\text{K}_{0}(\sqrt{s}\sqrt{r_{\text{m}k,p\text{D}}^{2}+ \alpha^{2}-2r_{\text{m}k,p\text{D}}\alpha\cos(\theta_{\text{f}k}-\theta_{\text {f}i})})\text{d}\alpha\right]+\\ &\frac{2\pi}{C_{\text{D}}}\left\{\sum\limits_{i=1}^{n}\overline{q }_{\text{f}k,i\text{D}}(s)\cdot[(v-i)\cdot\Delta L_{\text{fRD}}^{2}]+\frac{ \Delta L_{\text{fRD}}^{2}}{8}\cdot\overline{q}_{\text{f}k,p\text{D}}(s)-r_{ \text{m}k,p\text{D}}\sum\limits_{i=1}^{i=N_{k}}\overline{q}_{\text{f}k,i\text{ D}}(s)\Delta L_{\text{fRD}}\right\}=0,\end{split}  \tag{33}\]  
\[\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{N_{i}}\left[\overline{q}_{\text{fi,jD}} (s)\Delta L_{\text{fRD}}\right]=\frac{1}{s}, \tag{34}\]  
where \(1\leq k\leq n\), \(1\leq v\leq N_{k}\) and \(\Delta L_{\text{fID}}=L_{\text{fID}}/N_{i}\).