Context. The measured average velocities in solar and stellar spectral lines formed at transition region temperatures have been difficult to interpret. The dominant redshifts observed in the lower transition region naturally leads to the question of how the upper layers of the solar (and stellar) atmosphere can be maintained. Likewise, no ready explanation has been made for the average blueshifts often found in upper transition region lines. However, realistic three-dimensional radiation magnetohydrodynamics (3D rMHD) models of the solar atmosphere are able to reproduce the observed dominant line shifts and may thus hold the key to resolve these issues. Aims. These new 3D rMHD simulations aim to shed light on how mass flows between the chromosphere and corona and on how the coronal mass is maintained. These simulations give new insights into the coupling of various atmospheric layers and the origin of Doppler shifts in the solar transition region and corona. Methods. The passive tracer particles, so-called corks, allow the tracking of parcels of plasma over time and thus the study of changes in plasma temperature and velocity not only locally, but also in a co-moving frame. By following the trajectories of the corks, we can investigate mass and energy flows and understand the composition of the observed velocities. Results. Our findings show that most of the transition region mass is cooling. The preponderance of transition region redshifts in the model can be explained by the higher percentage of downflowing mass in the lower and middle transition region. The average upflows in the upper transition region can be explained by a combination of both stronger upflows than downflows and a higher percentage of upflowing mass. The most common combination at lower and middle transition region temperatures are corks that are cooling and traveling downward. For these corks, a strong correlation between the pressure gradient along the magnetic field line and the velocity along the magnetic field line has been observed, indicating a formation mechanism that is related to downward propagating pressure disturbances. Corks at upper transition region temperatures are subject to a rather slow and highly variable but continuous heating process. Conclusions. Corks are shown to be an essential tool in 3D rMHD models in order to study mass and energy flows. We have shown that most transition region plasma is cooling after having been heated slowly to upper transition region temperatures several minutes before. Downward propagating pressure disturbances are identified as one of the main mechanisms responsible for the observed redshifts at transition region temperatures.