The overall performance of a tokamak strongly depends on phenomena that take place in a thin region between the main plasma and the vessel wall, which is denoted as tokamak boundary. In fact, the formation of transport barriers in this region can significantly improve plasma confinement and, therefore, the tokamak fusion performance. In addition, the tokamak boundary determines the peak heat flux to the wall, an essential quantity for the design and the operation of fusion power plants, as well as the level of impurities in the core, the removal of fusion ash and the dynamics of neutral particles. The dynamics in the plasma boundary is strongly nonlinear and characterized by a wide range of length and time scales as well as by a complex magnetic field geometry that may feature one or more nulls of the poloidal magnetic field. Large-scale, three-dimensional turbulence simulations are therefore often required to disentangle the complex physical mechanisms that govern this region. The thesis is focused on the analysis of the different turbulent transport regimes present in the plasma boundary as they appear from three-dimensional, flux-driven, global, two-fluid turbulence simulations carried out by using the GBS code, which is significantly extended here to allow the self-consistent simulation of the plasma dynamics coupled to a kinetic single-species neutral model in arbitrarily complex magnetic geometries. Considering single-null magnetic configurations, three turbulent transport regimes are identified: (i) a regime of suppressed turbulent transport at low values of collisionality and large values of heat source, (ii) a regime of developed turbulent transport at intermediate values of collisionality and heat source, and (iii) a regime of very large turbulent transport at high value of collisionality and density, which can be associated to the crossing of the density limit. By leveraging the results of GBS simulations, theory-based scaling laws of the pressure and density decay lengths in the near and far scrape-off layer are derived in the developed transport regime from a balance among heat source, turbulent transport across the separatrix and parallel losses at the vessel wall. The theoretical scaling of the pressure decay length in the near scrape-off layer is successfully validated against a multi-machine database of SOL width measurements at the outer target. By carefully analysing the transition to the regime of large turbulent transport, we show that the density limit can be explained by an enhancement of turbulent transport at the tokamak boundary when the density increases. This analysis leads to a theory-based scaling law of the maximum edge density achievable in tokamaks, which is in better agreement with a multi-machine database than the widely used Greenwald empirical scaling, thus significantly improving our understanding and predictive capability of the density limit, with important implications for the design and the operation of future fusion power plants. The thesis concludes by presenting the first turbulent simulations carried out in various snowflake magnetic configurations, which are used to investigate the effect of turbulence and equilibrium flow on the heat flux distribution among the different strike points.