This article compares the accuracy of return value estimates from stationary and non-stationary extreme value models when the data exhibits covariate dependence. The non-stationary covariate representation used is a penalised piecewise-constant (PPC) model, in which the data are partitioned into bins defined by covariates and the extreme value distribution is assumed to be homogeneous within each bin. A generalised Pareto model is assumed, where the scale parameter can vary between bins but is penalised for the variance across bins, and the shape parameter is assumed constant over all covariate bins. The number and sizes of covariate bins must be defined by the user based on physical considerations. Numerical simulations are conducted to compare the performance of stationary and non-stationary models for various case studies, in terms of quality of estimation of the -year return value over the full covariate domain. It is shown that a non-stationary model can give improved estimates of return values, provided that model assumptions are consistent with the data. When the data exhibits non-stationarity in the generalised Pareto tail shape, the use of non-stationary model assuming a constant shape parameter can produce biases in return values. In such cases, a stationary model can give a more accurate estimate of return value over the full covariate domain as only the most extreme observations (regardless of covariate) are used to estimate tail shape. In other cases, the assumption of a stationary model will ignore key features of the data and be less reliable than a non-stationary model. For example, if a relatively benign covariate interval exhibits a long (or heavy) tail, extreme values from this interval may influence the T-year return value for very large T. However the sample of peaks over threshold, with high threshold, used to estimate a stationary model in this case may not include sufficient observations from this interval to estimate the return value adequately.