Large eddy simulations of the flow around a square-back road vehicle (e.g. a prototype of a small van) have been performed. Three different configurations are considered. The baseline geometry has a regular square back with sharp edges. The second configuration has an appendix on the tail consisting of a solid protrusion of 50 mm and a tapered angle of 9 deg. The third model has four plates that form a cavity attached to the base. The numerical technique is based on the immersed boundary approach, which allows the use of underlying Cartesian grids for the analysis of complex configurations by prescribing suitable body forces to enforce the no-slip boundary condition on solid surfaces. Grids with up to 10 million cells were used to study the effect of the Reynolds number Re and of the subgrid-scale model. Comparisons with the experimental data show that time-averaged quantities, as well as flow dynamics, are accurately predicted. In particular, mean velocity profiles, total drag, base pressure, low-frequency axial wake pumping, and high- frequency shear layer instabilities are very close to the measured values. At Re = 2 x 10⁴, the most evident feature of the flow around the baseline configuration is the formation of a large separation bubble at the base of the body, which is responsible for most of the drag. This region is confined by four shear layers generated at the edges of the base, which force the recirculation region. The numerical results show a secondary recirculation region on the ground in the rear vehicle region, which vanishes at higher Reynolds numbers. The predicted drag reduction efficiency of the second and the third configuration is also in very good agreement with the experiments. The main difference to the baseline configuration is the absence of the ground recirculation at low Reynolds number. Time histories of lift and side-force coefficients show that the plates strongly reduce the flow unsteadiness with respect to the other configurations. The boat tail geometry produces the most drag reduction, due to the pressure recovery and the smaller recirculation region. A comparison of the results produced with two different subgrid scale stress models, the Smagorinsky model and a dynamic model, shows, that the dynamic model is clearly superior because it allows an adaptive distribution of turbulent viscosity, which is especially important in complex geometry flows.