Both hydraulic fracture permeability and matrix permeability play important roles in the life of gas production. However, the matrix permeability has not yet been understood fully because gas flow within the matrix undergoes a transition from viscous flow to slip flow and Knudsen diffusion. Traditional Darcy law cannot be utilized directly to represent this phenomenon. Thus, understanding the gas flow inside the matrix and how the matrix permeability evolves during gas depletion has been a major research challenge. In this study, an apparent permeability model for shale matrix is developed to reflect the combined impact of flow regimes and effective stress. Flow regimes are defined by the Knudsen number while the effect of effective stress is expressed by the variable porosity and intrinsic permeability based on poroelasticity theory. This model is verified against the experimental data of the permeability for shale plug samples. Then, the apparent permeability is used to couple shale deformation and gas flow within the matrix. The fully coupled model is implemented and solved by Comsol Multiphysics based on finite element method. Numerical simulation results demonstrate that both flow regimes and effective stress have significant impacts on the apparent permeability of the matrix. The apparent permeability for shale matrix undergoes a slight increase at the early time of production and becomes bigger at the later time because of slip flow and Knudsen diffusion. The slippage effect has a critical impact on the apparent permeability when gas pressure drops to a low magnitude and the Knudsen number would be higher than 0.1. However, the intrinsic permeability experiences a declining trend as the effective stress increases during gas production. Sensitivity analyses indicate that both initial intrinsic permeability and initial porosity have more significant effects on the apparent permeability than the bulk modulus of shale matrix.