We investigate a connection between spatial statistics and statistical physics to obtain new covariance functions with direct physical interpretation for spatial random fields. These covariance functions are based on the exponential Boltzmann-Gibbs representation and use an energy functional to represent interactions between the values of the random field at different points in space. This formulation results in closed-form generalized covariance functions, which display infinite variance in Euclidean spaces of dimension larger than one. We propose regularization schemes in real and reciprocal (spectral) space that lead to well-behaved covariance structures. The real-space regularization parameter allows a continuous interpolation between the Boltzmann-Gibbs covariance and the exponential covariance. We also propose discretized approximations on regular grids, and we show that they represent reparametrized versions of the well-known Besag and Leroux lattice models. We then discuss parameter estimation and spatial prediction for the regularized Boltzmann-Gibbs covariance model in two dimensions. We recommend using the pairwise difference likelihood that combines satisfactory estimation performance and good scalability with many observation points. The predictive performance of the regularized covariance function is assessed by means of cross-validation statistics. Irregularlyspaced samples from the Walker Lake dataset are used, and spatial prediction is conducted by means of ordinary kriging. The regularized Boltzmann-Gibbs covariance yields improved predictive performance compared to the exponential covariance model.