Enhanced gas recovery (EGR) is believed to be a promising technology to improve the production of shale gas reservoirs and simultaneously reduce the emissions of greenhouse gas via the injection (sequestration) of carbon dioxide, to which great effort has been devoted by scholars. However, traditional investigations are generally limited to the ideal model of nanochannels and statistic characterization of competitive adsorption, neglecting the nanoporous structure of the kerogen matrix and the complex dynamic behavior during the EGR process. In this work, we present a comprehensive study of the EGR process in a realistic kerogen pore network (matrix) which is obtained from the artificial pyrolysis of bulk kerogen through reactive force field molecular dynamics (ReaxFF MD) simulations. The influence of pore properties (e.g., porosity) of the kerogen matrix under different maturities, and the proportions (i.e., methane and carbon dioxide) of injection gas with various injection pressures are revealed and meticulously discussed. In addition, the underlying mechanisms including diffusion and displacement effects behind the EGR process are analyzed by combining them with with particle trajectory capture technology. In particular, based on the MD simulation results, an analytical model to depict the dynamic recovery process in the kerogen matrix is proposed by coupling consideration of recovery time and capacity, which are examined against the simulation and experimental data. The hybrid recovery strategy is developed by utilizing the advantages of depressurization and gas-injection recoveries to achieve the optimization of both recovery time and capacity. The insights acquired from this work would be helpful for efficient exploitation of shale gas reservoirs and pave the way to capture the realistic EGR processes within the kerogen matrix from molecular and theoretical perspectives.