Natural convection gaseous slip flows in open-ended vertical parallel-plate microchannels with symmetric wall heat fluxes are numerically investigated. A second-order model, including thermal creep effects, is considered for velocity slip and temperature jump boundary conditions with variable thermophysical properties. Simulations are performed for wide range of Rayleigh numbers from 5 x 10(- 6) to 5 x 10(- 3) in the continuum to slip flow regime. The developing and fully developed solutions are examined by solving the Navier-Stokes and energy equations using a control volume technique. It is found that the second-order effects reduce the temperature jump and the slip velocity, whereas thermal creep strongly increases the slip velocity in both developing and fully developed regions. Moreover, the rarefaction effects increase the flow and heat transfer rates considerably, while decreasing the maximum gas temperature and friction coefficient as compared to the continuum limit. It was also shown that the axial temperature variations of the gas layer adjacent to the wall in the modeling of the thermal creep are of paramount importance and neglecting these variations, which is common in literature, leads to unphysical velocity and temperature distributions.