We propose a simple piecewise model for a sample of peaks-over-threshold, nonstationary with respect to multidimensional covariates, and estimate it using a carefully designed and computationally efficient Bayesian inference. Model parameters are themselves parameterized as functions of covariates using penalized B-spline representations. This allows detailed characterization of non-stationarity extreme environments. The approach gives similar inferences to a comparable frequentist penalized maximum likelihood method, but is computationally considerably more efficient and allows a more complete characterization of uncertainty in a single modelling step. We use the model to quantify the joint directional and seasonal variation of storm peak significant wave height at a northern North Sea location and estimate predictive directional–seasonal return value distributions necessary for the design and reliability assessment of marine and coastal structures.