## 1 Introduction  
Shale gas has emerged as a promising resource to fulfill global energy demand, primarily attributed to notable discoveries of reserves and the combined application of horizontal drilling techniques and multi-stage hydraulic fracturing methods [1, 2]. According to data from the U.S. Energy Information Administration (EIA), the United States possesses technically recoverable shale gas resources amounting to 862 trillion cubic feet [3]. These resources are distributed across diverse shale formations, including the Marcellus, Utica, Barnett, and Haynesville. These shale formations are recognized for their low porosity and ultra-low permeability. Nevertheless, advancements in horizontal drilling and hydraulic fracturing technologies have rendered the extraction of natural gas from these formations economically feasible [4, 5]. This entails the initial vertical drilling of a well, followed by its subsequent horizontal orientation within the shale formation. Sequential simulations of multiple stages are then implemented along the horizontal wellbore through the injection of high-pressure fluids, typically comprising water, sand, and additives. This process leads to the creation of a stimulated reservoir volume, thereby amplifying the contact areas between the wellbore and the shale formation. Consequently, the horizontal wells with multistage hydraulic fracturing enhance gas drainage efficiency and augment overall gas production. In the year 2022, dry natural gas production from the shale formations in the United States was estimated at 28.5 trillion cubic feet due to the combination of horizontal drilling with multistage hydraulic fracturing, equating to roughly 80% of the total dry natural gas production [6]. To enhance shale gas recovery from hydraulically fractured reservoirs, the prediction of production at a specified time is imperative for selecting the optimal completion method and well operational parameters, as well as for refining gas field development strategies [7; 8; 9]. Conventional methods for production forecasting, such as decline curve analysis (DCA) and numerical reservoir simulation, are extensively employed in the oil and gas sector to foresee well production performances from hydrocarbon reservoirs. DCA, a prevalent empirical approach, relies on historical production data to anticipate future production outcomes [10; 11]. This method assumes that the production decline of a well or a reservoir follows a certain mathematical function. Through fitting the decline curve to the available production data, essential parameters like the initial production rate, decline rate, and ultimate recovery can be estimated. DCA offers a straightforward and expedient estimation of forthcoming production; however, it operates under the assumption of boundary-dominated flow and may not accurately capture intricate reservoir behavior or shifts in production mechanisms over time. Conversely, numerical reservoir simulation considers reservoir heterogeneities, well configurations, fluid properties, and other factors exerting influence on hydrocarbon production [12; 13; 14]. This approach permits the integration of intricate reservoir dynamics, encompassing fluid flow dynamics, pressure variations, and reservoir geometry. For instance, the reservoir simulations adeptly incorporate various intricate mechanisms governing shale gas transport in the production prediction, including Fick diffusion and non-Darcy flow, alongside the geometry of the simulated reservoir volume [15; 16]. Nonetheless, the reservoir simulations require substantial data input encompassing reservoir properties, well specifics, and fluid characteristics. The accuracy of production predictions heavily relies on the quality of the input data and the reliability of the reservoir models. Moreover, the simulations usually demand considerable computational resources and entail time-intensive model construction and execution processes.