Consider a symmetric function (Mercer Kernel) $K : T \times T \rightarrow \mathbb{R}$ and define an operator $H_k: L^2(T,\nu) \rightarrow L^2(T,\nu)$ where $H_kf(x) = \int_X K(x,y)f(y)d\nu(y)$.
Further suppose we have $$ f(x) = \sum_i \sqrt{\lambda_i}Z_i\phi_i(x) $$ with $\phi_1,\phi_2,\ldots$ a countably infinite orthonormal basis of eigenfunctions for the space spanned by $H_k$ (with associated $\lambda_i$, we know all this from Mercer's theorem).
Next consider $V_k$ to be the space spanned by the first $k$ eigenfunctions $\{\phi_1,\ldots,\phi_k\}$ of $H_k$. As $f(x)$ is defined to be the sum of elements in the (possibly infinite) orthonormal basis, we can define the projection of $f$ onto the space $V_k$ as $$ f_k = \sum_{i=1}^k \sqrt{\lambda_i}Z_i\phi_i(x) $$
I believe $f_k$ is a finite-dimensional normal RV, as it is a sum of independent normal random variables. But what is the variance of $f_k$?