AbstractFor a linear second order elliptic partial differential operator A : V → V′, we consider the boundary value problems Au = f with stationary Gaussian random data f over the dual V′ of the separable Hilbert space V in which the solution u is sought. The operator A is assumed to be deterministic and bijective. The unique solution u = A−1f is a Gaussian random field over V. It is characterized by its mean field Eu and its covariance Cu ∈ V ⊗ V. For a class of Gaussian data f with piecewise analytic covariance kernels Cf ∈ V′ ⊗ V′, we characterize analytic regularity of the covariance Cu of the Gaussian solution u in families of countably normed, weighted Sobolev spaces. To this end, we investigate shift theorems for the (nonhypoelliptic) deterministic tensor PDEs (A ⊗ A)Cu = Cf proposed by Schwab and Todor, (Numer Math 95 (2003), 707–734) for the computation of covariance Cu of the random solution u. The nonhypoelliptic nature of A ⊗ A implies that sing supp(Cu) is in general strictly larger than sing supp(Cf). For a model problem, analyticity and singular support of the solution is characterized completely. Based on our regularity results, we outline an hp‐finite element strategy (Pentenrieder, Ph.D thesis, ETHZurich, Dissertation No. 18729, 2009; Pentenrieder and Schwab, Research Report 2010‐08, Seminar for Applied Mathematics, ETH Zürich, Submitted) to approximate Cu stemming from covariances of stationary Gaussian data f. In the second part (Pentenrieder and Schwab, Research Report 2010–08, Seminar for Applied Mathematics, ETH Zürich, Submitted) of this work, we prove that this discretization gives exponential rates of convergence of the FE approximations, in terms of the number of degrees of freedom. © 2011 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq, 2012