Molecular diffusion represents one of the critical mass transfer mechanisms of water-gas-oil in subsurface shale porous media with complex fracture networks (CFNs). The diffusion coefficient is a crucial parameter governing the mass transfer by molecular diffusion. Nonetheless, there remains insufficient evidence to ascertain the optimal empirical correlation for predicting diffusion coefficients. This paper, therefore, presents an improved Embedded Discrete Fracture Model (EDFM) to determine the optimal empirical correlation that calculates the diffusion coefficients. The improved EDFM firstly integrated four widely used empirical correlations. Taking the experimental results as a reliable reference, the optimal empirical correlation was identified by comparing the calculation predicted by the improved EDFM. The comparison results reveal that the extended Sigmund empirical correlation is the optimal method for calculating diffusion coefficients in subsurface shale porous media with CFNs. Based on that, an in-depth investigation was performed into the mass transfer of water-gas-oil considering molecular diffusion during hydrocarbon gas huff ’n’ puff. The results indicate that molecular diffusion reduces the convective mass transfer rate (C-MTR) of injected components from connected natural fractures (CNFs) to the matrix during both the injection and soaking stages. However, the total mass transfer rate (T-MTR) is enhanced. For non-injected components, molecular diffusion has a minimal effect on the C-MTR. However, during the production stage, molecular diffusion enhances the T-MTR of lighter non-injected components from the matrix to the CNFs. This paper establishes a theoretical foundation for comprehensive investigation of mass transfer of water-gas-oil considering molecular diffusion in subsurface shale porous media with CFNs.