Let $ψ:{\mathcal{D}}\rightarrow{\mathbf{R}}$ be a harmonic function such that $Δψ(x)=0$ for all $x\in\mathcal{D}\subset{\mathbf{R}}^{n}$. There are then many well-established classical results:the Dirichlet problem and Poisson formula, Harnack inequality, the Maximum Principle, the Mean Value Property etc. Here, a 'noisy' or random domain is one for which there also exists a classical scalar Gaussian random field (GRF) ${\mathscr{J}(x)}$ defined for all $x\in{\mathcal{D}}$ or $x\in\partial {\mathcal{D}}$ with respect to a probability space $[Ω,\mathcal{F},{\mathrm{I\!P}}]$. The GRF has vanishing mean value $\mathbf{E}[\![\mathscr{J}(x)]\!] = 0$ and a regulated covariance ${\mathbf{E}}[\![{\mathscr{J}(x)} \otimes {\mathscr{J}(y)}]\!] = αJ(x,y;ξ)$ for all $(x,y)\in{\mathcal{D}}$ and/or $(x,y)\in{\partial\mathcal{D}}$, with correlation length $ξ$ and ${\mathbf{E}}[\![{\mathscr{J}(x)} \otimes {\mathscr{J}}(x)]\!] = α<\infty$. The gradient $\nabla{\mathscr{J}(x)}$ and integral $\int_{\mathcal{D}}{\mathscr{J}}(x) dμ(x)$ also exist on ${\mathcal{D}}\bigcup\partial\mathcal{D}$. Harmonic functions and potentials can become randomly perturbed GRFs of the form $\overline{ψ(x)}=ψ(x)+λ{\mathscr{J}}(x)$. Physically, this scenario arises from noisy sources or random fluctuations in mass/charge density, noisy or random boundary/surface data; and introducing turbulence/randomness into smooth fluid flows, steady state diffusions or heat flow. This leads to stochastic modifications of classical theorems for randomly perturbed harmonic functions and Riesz and Newtonian potentials; and to stability estimates and bounds for the growth and decay of their volatility and moments.