$\newcommand{\f}{\mathbf f}\newcommand{\e}{\varepsilon}\newcommand{\y}{\mathbf y}\newcommand{\0}{\mathbf 0}$If $f \sim \mathcal {GP}(0, k)$ and we have sample points $x_1,\dots,x_n$ then letting $\f = (f(x_1), \dots, f(x_n))^T$ we have
$$
\f \sim \mathcal N(\0, K)
$$
with $K_{ij} = k(x_i,x_j)$ as the kernel matrix. If we observe
$$
y_i = f(x_i) + \e_i
$$
so
$$
\y = \f + \e
$$
where $\e \sim \mathcal N(\0, \sigma^2_\e I) \perp f$ then
$$
\y \sim \mathcal N(\0, K + \sigma^2_\e I).
$$
Now for a new point $x_s$ we have two options. If we are just interested in the posterior distribution of a point of the Gaussian process itself, i.e. $f_s := f(x_s)$, we'll have
$$
{\f \choose f_s} \sim \mathcal N\left({\0 \choose 0}, \begin{bmatrix}K & k_s^T \\ k_s & k_{ss}\end{bmatrix}\right)
$$
so
$$
{\y \choose f_s} = {\f \choose f_s} + {\e \choose 0} \sim \mathcal N\left({\0 \choose 0}, \begin{bmatrix}K + \sigma^2_\e I & k_s^T \\ k_s & k_{ss}\end{bmatrix} \right)
$$
which means the posterior variance of $f_s\mid \y$ is
$$
\text{Var}(f_s\mid \y) = k_{ss} - k_s^T(K + \sigma^2_\e I)^{-1}k_s.
$$
If however we are interested in the posterior distribution of a new noisy data point $y_s = f_s + \e_s$ then we'll have
$$
{\y \choose y_s} = {\f \choose f_s} + {\e \choose \e_s} \sim \mathcal N\left({\0 \choose 0}, \begin{bmatrix}K & k_s^T \\ k_s & k_{ss}\end{bmatrix} + \sigma^2_\e I_{n+1}\right).
$$
so posterior variance of $y_s\mid \y$ is
$$
\text{Var}(y_s\mid \y) = (k_{ss} + \sigma^2_s) - k_s^T(K + \sigma^2_\e I)^{-1}k_s.
$$
Thus it depends on whether or not we're talking about the posterior variance of the GP itself versus new data points coming from the noisy process. If you want to get a sense of the variation of the GP $f$ then you can look at $f_s \mid \y$ but if you're interested in what noisy data will look like coming from it then $y_s\mid \y$ will be more interesting.
In linear regression this is analogous to looking at the estimated mean for a new point $\text{Var}(x_s^T\hat\beta)$ versus looking at the variance of a new data point noise and all $\text{Var}(x_s^T\hat\beta + \e_s)$. Both can be reasonable but it depends what we're trying to do. Often though in my experience we care more about the mean, i.e. $f$ or $\hat\beta$, so usually that's what I go with, rather than the predictive distribution for a new data point.
Here's an example showing the posterior mean with 1.96 times the standard deviation of the variance for just the GP itself vs for new data. The underlying GP is well captured by the red line while the data is well described by the blue line, which just has the extra $\sigma^2_\e$ term.
