---

Spatial regression models are used to analyze quantitative spatial responses based on linear and non-linear relationships with explanatory variables. Usually, the spatial correlation of responses is modeled with a Gaussian random field based on a multivariate normal distribution. However, in practice, we encounter skewed responses, which are analyzed using skew-normal distributions. Closed skew-normal distribution is one of the extended families of skew-normal distributions, which has similar properties to normal distributions. This article presents a hierarchical Bayesian analysis based on a flexible subclass of closed skew-normal distributions. Given the time-consuming nature of Monte Carlo methods in hierarchical Bayes analysis, we have opted to use the variational Bayes approach to approximate the posterior distribution. This decision was made to expedite the analysis process without compromising the accuracy of our results. Then, the proposed model is implemented and analyzed based on the real earthquake data of Iran.