The grand-canonical Monte Carlo (GCMC) and molecular dynamics (MD) simulations are performed to investigate the methane aqueous fluids in Na-montmorillonite clay interlayer under near-surface geological temperature and pressure conditions (T = 300 K and P = 20-50 bar). The chemical potentials of water and methane under these T/P conditions are calculated using the Widom's insertion method. These chemical potentials are used in the GCMC simulations to determine the contents of different species in the clay interlayer, especially in those that correspond to the equilibrium stable spacing distances. Simulation results show that initial clay swelling is dominated by water adsorption into the clay interlayer, followed by the intercalation of methane as the basal spacing increases. However, it is found that this methane intercalation process is strongly influenced by the relative humidity and the total gas pressure of the system. High relative humidity may facilitate water molecules entering the clay interlayer region and inhibit the intercalation of methane molecules. MD simulations show that sodium ions are fully hydrated by water molecules and clay surface oxygen atoms, while methane molecules are not fully coordinated. This situation is attributed to the less water content in clay interlayer and the subsequent formation of methane dimer or trimer clusters due to the hydrophobic nature of small hydrocarbon molecules.