In this paper, natural convection melting in a square cavity with gradient porous media is numerically studied at pore-scale level by adopting the lattice Boltzmann method. To generate the gradient porous media, a Monte Carlo technique based on the random sampling principle is used. The effects of several factors, such as Rayleigh number, gradient porosity structure, gradient direction, and particle diameters on natural convection melting are investigated in detail. Based on the numerical data, it is observed that the thermal performance of the gradient porous media always depends on the Rayleigh number and, specifically, as the Rayleigh number is set to 106, the total melting time obtained for the case of the negative gradient porous media is always shorter than the cases of positive gradient and uniform porous media. However, if the Rayleigh number is equal to 104, at which the heat transfer is dominated by the heat conduction, it is noted that the performance of the positive gradient porous media is better than the other cases. To have a better understand on this point, various simulations are also performed and we found that there usually exists a critical value of Rayleigh number to determine the thermal performance of the gradient porous media. Moreover, our numerical results also show that the influence of the particle diameter on the liquid fraction is insignificant as Rayleigh number is set to 104, while it has a great impact on the liquid fraction when Rayleigh number equals 106.