Fitting spatial models with a Gaussian random field as spatial random effect poses computational challenges for Markov Chain Monte Carlo (MCMC) methods, primarily due to two factors: computational speed and convergence of chains for the hyperparameters. To deal with this, a Gaussian random field can be approximated by a Gaussian Markov random field using stochastic partial differential equations. This methodology is commonly used in “latent Gaussian models”, where the inference is done by the Integrated Nested Laplace Approximations, but rarely used in an MCMC method. In this contribution, we evaluated different parameterizations of the approximated Gaussian random field, specifically using the Hamiltonian Monte Carlo algorithm of the Stan software. A simulation study demonstrated that models using the hyperparameters ρ and σu were better able to estimate the values used to simulate the spatial random field. Their speed computation were faster compared to models parameterized with κ and τ. In real data application, the index of relative abundance estimated for Pollock indicates similar trends for the six models proposed. However, models incorporating ρ and σu demonstrated faster computation compared to those utilizing κ and τ, corroborating the results found in the simulation. Even more important, none of these models encountered convergence issues, as indicated by the Rhat statistic.