A novel first-principles self-consistent model that couples plasma and neutral atom physics suitable for the simulation of turbulent plasma behaviour in the tokamak edge region has been developed [1] and implemented in the GBS code [2]. While the plasma is modelled by the drift-reduced two fluid Braginskii equations, a kinetic model is used for the neutrals, valid in short and in long mean free path scenarios. The model includes ionization, charge-exchange, recombination, and elastic collisional processes. The model was used to study the transition from the sheath to the conduction limited regime by increasing the plasma density in the system. We compared the simulation results with the predictions of an expanded and refined two-point model, estimating the drop of electron and ion temperature along the magnetic field lines in the SOL. The model is now being applied to investigate the impact of neutrals on turbulent edge features, such as the broad shoulder observed in the far SOL at high plasma density. References [1] C. Wersal and P. Ricci 2015 Nucl. Fusion 55 123014 [2] P Ricci et al 2012 Plasma Phys. Control. Fusion 54 124047