Given a random scalar field $u(x,y)$ for $x,y \in D \subset \mathbb{R}^2$ where $D$ is a rectangular box centered at $(0,0)$, the value of $u(x,y)$ at an arbitrary point $(x,y)$ is uniformly distributed in $[0,1]$, i.e. similar to a white noise but is uniformly distributed rather than normally distributed. The field is convolved with a Gaussian function (a deterministic function, not a random distribution) in the form of $g(x,y)=a \exp(-(x^2+y^2)/2\sigma^2)$ to give another random scalar field $f(x,y)$, i.e. $$ f(x,y) = u(x,y) * g(x,y), $$ or in terms of Fourier transforms $$ f(x,y) = \mathcal{FT}^{-1}[\mathcal{FT}(u)\mathcal{FT}(g)](x,y). $$ If the random scalar field $u(x,y)$ is sampled multiple times to obtain $u_i(x,y)$ for $i=1,2,...,N$, such that we obtain different random scalar fields $f_i(x,y)$, what can we say about the statistical properties (e.g. statistical moments) of an individual field $f_i(x,y)$ and the average field $f_{mean}(x,y):=\frac{1}{N}\sum_{i=1}^N f_i(x,y)$ for $N \rightarrow \infty$? Are the autocorrelations $\langle f_{i}(\vec{x_0})f_{i}(\vec{x_0}+\vec{x})\rangle$ and $\langle f_{mean}(\vec{x_0})f_{mean}(\vec{x_0}+\vec{x})\rangle$ in Gaussian form? Can we ensure that $f(x,y) > 0$, at least in most of the cases?
The primary purpose here is to generate a non-negative correlated random scalar field, such that the autocorrelation of the strengths of the field between two points is related by a Gaussian function (i.e. the closer the two points are, the more correlated the field strengths at these two points are; the standard deviation of the Gaussian is a typical length scale characterizing the smoothness of the random scalar field). Since non-negativity is preferred, uniform distribution is used here, in the hope that the autocorrelation is not too non-Gaussian.