Direct survey estimates for small areas are likely to yield unacceptably large standard errors due to the small sample sizes in the areas. This makes it necessary to use models to“borrow strength" from related areas to find more reliable estimate for a given area or, simultaneously, for several areas. For instance, in many applications, data on related multiple characteristics and auxiliary variables are available. Thus, multivariate modeling of related characteristics with multiple regression can be implemented. However, while model-based small area estimates are very useful, one potential difficulty with such estimates when models are used is that the combined estimate from all small areas does not usually match the value of the single estimate on the large area. Benchmarking is done by applying a constraint to ensure that the“total" of the small areas matches the“grand total". Benchmarking can help to prevent model failure, an important issue in small area estimation. It can also lead to improved prediction for most areas because of the information incorporated in the sample space due to the additional constraint. We describe both the univariate and multivariate Bayesian nested error regression models and develop a Bayesian predictive inference with a benchmarking constraint to estimate the finite population means of small areas. Our models are unique in the sense that our benchmarking constraint involves unit-level sampling weights and the prior distribution for the covariance of the area effects follows a specific structure. We use Markov chain Monte Carlo procedures to fit our models. Specifically, we use Gibbs sampling to fit the multivariate model; our univariate benchmarking only needs random samples. We use two datasets, namely the crop data (corn and soybeans) from the LANDSAT and Enumerative survey and the NHANES III data (body mass index and bone mineral density), to illustrate our results. We also conduct a simulation study to assess frequentist properties of our models.