A model for the simulation of the tokamak boundary is described, which is mass-conserving to leading order in rho s / R 0 (where rho s is the ion sound Larmor radius and R-0 is the tokamak major radius at the magnetic axis). The model integrates a three-dimensional plasma turbulence model, based on the two-fluid drift-reduced Braginskii equations, a neutral model, based on the discretization of the Boltzmann equation integrated along the neutral path, and proper boundary conditions at the vessel walls. The implementation of the model in the GBS code [Ricci et al., Plasma Phys. Control. Fusion 54, 124047 (2012)] is described and a technical proof of mass conservation is obtained. The simulation results highlight the impact of mass conservation on the plasma and neutral profiles, and the formation of the plasma density profile is explained in terms of plasma and neutral particle fluxes.