Statistical post-processing is an indispensable tool for providing accurate weather forecasts and early warnings for weather extremes. Most statistical post-processing is univariate, with dependencies introduced via use of an empirical copula. Standard empirical methods take a dependence template from either the raw ensemble output (ensemble copula coupling, ECC) or the observations (Schaake Shuffle, SSh). There are drawbacks to both methods. In ECC it is assumed that the raw ensemble simulates the dependence well, which is not always the case (e.g. 2-meter temperature in The Netherlands). The Schaake Shuffle is not able to capture flow dependent changes to the dependence and the choice of observations is key. Here we compare a reshuffled standard ensemble model output statistics (EMOS) approach with two multivariate bias adjustment approaches that have not been used before in a post-processing context: 1) the multivariate bias correction with N-dimensional probability density function transform (MBCn) and 2) multivariate ranks that are defined with optimal assignment (OA). These methods have the advantage that they are able to explicitly capture the dependence structure that is present in the observations. We apply ECC, the Schaake Shuffle, MBCn and OA to 2-m and dew point temperature forecasts at seven stations in The Netherlands. Forecasts are verified with both univariate and multivariate methods, and using a heat index derived from both variables, the wet-bulb globe temperature (WBGT). Our results demonstrate that the spatial and inter-variable dependence structure is more realistic in MBCn and OA compared to ECC or the Schaake Shuffle. The variogram score shows that while ECC is most skilful for the first two days, at moderate lead times MBCn is most skilful and at the longest lead times OA is more skilful than both ECC and MBCn. Overall, we highlight the importance of considering the dependence between variables and locations in the statistical post-processing of weather forecasts.