After fracturing operation, hydraulic fractures and induced fractures are created within the shale reservoir. A lot of treatment water is stored in fractures network and flow back into the surface during the gas recovery process. The gas production performance is affected by the water flowback because two phase flow occurs within fractures zone. For the created reservoir scale, we propose a multiscale- multiphase simulation model, which defines the whole domain as three sections. Section A contains the organic and inorganic matrix, which stores both the free gas and adsorbed gas. Flow processes are defined in the components of inorganic minerals and kerogens, respectively. For the section B and C, gas phase and water phase are existed together. Under this framework, a set of partial differential equations are derived to define various liquid transport processes: (1) gas flow in the kerogen system of matrix; (2) gas flow in the inorganic system of matrix; (3) gas-water two phase flow in fractures zone and (4) gas-water two phase flow in the hydraulic fracture system. Dynamic permeability models and mass exchanges between them are coupled for all systems. The model was verified against field production data from the Barnett Shale. Model simulation results show that flowback of treatment water can significantly affect the gas production rate at the early stage. Firstly, the increase of maximum water relative permeability can raise the water flowback rate and gas production rate but increasing non-wetting phase entry pressure will decrease the fluids flow rate. Secondly, the impact of fractures zone width on gas production performance is unstable and increasing initial water saturation can increase the water flowback rate but decrease gas production rate. Overall, the dynamic performances of water phase within fractures zone have significant impact on the short and long time shale gas recovery.