The sample mean of a multivariate normal vector $\mathbf{X}=(X_1, X_2, \ldots, X_n)$ is a function of
$$M = X_1+X_2+\cdots X_n$$
and the sample variance is a function of the residual vector with components
$$Z_i = -X_1 - X_2 - \cdots - X_{i-1} + (n-1)X_i - X_{i+1} - \cdots - X_n,$$
$i=1, 2, \ldots, n$.
Let $\Sigma$ be the covariance matrix of $\mathbf{X}$. Write $\sigma_i$ for the sum of column (or row) $i$ of $\Sigma$, $\sigma_i = \Sigma_{1i} + \Sigma_{2i} + \cdots + \Sigma_{ni}$, and let $\sigma$ be the sum of all the entries of $\Sigma$. We may compute
$$\operatorname{Cov}(M, Z_i) = n\sigma_i - \sigma.$$
Because both $M$ and $Z_i$ are linear combinations of multivariate Normal variables, they are jointly Normal, whence they are independent if and only if their covariance is zero. Consequently $M$ is independent of all the $Z_i$ if and only if
$$n\sigma_1 = n\sigma_2 = \cdots = n\sigma_n = \sigma.$$
In other words, equality of the column sums guarantees independence of the mean and the components of the sample variance, whence it will guarantee independence of the mean and the sample variance itself.
Although the converse is not true--it is possible for $M$ not to be independent of the $Z_i$, yet for $M$ to be independent of the sample mean--this requires exceptional circumstances. In almost all cases, inequality of the column sums creates a dependence between the sample mean and sample standard deviation.
By definition, in a stationary process the covariances $\Sigma_{ij}$ may depend only on $i-j$. Although this does not guarantee the column sums are all equal, for large $n$ and a covariance that decays sufficiently rapidly with $|i-j|$, it will approximately be true, because in the limit the column sums are all equal:
$$\sigma_i = \sum_{j=-\infty}^\infty \Sigma_{ji} =\sum_{j=-\infty}^\infty \Sigma_{jk} =\sigma_k.$$
All that is required is the convergence of these sums.
A good way to see the dependence in the scatterplot is to render the points more carefully. When they are made semitransparent, you can see the underlying density better. A lowess smooth helps demonstrate a variation in the standard deviation with the mean in this example where $n=8$ and the column sums of $\Sigma$ vary appreciably.

Here is the R
code that generated it.
library(MASS) # mvrnorm()
set.seed(17)
n <- 5e4 # Simulation size
d <- 8 # Dimension
k <- 4 # Size of upper block of Sigma
rho <- 0.99 # Correlation in upper block
mu <- rep(0, d) # Mean
Sigma <- outer(1:d, 1:d, function(i,j) ifelse(i <= k & j <= k, rho^abs(i-j), i==j))
colSums(Sigma)
x <- mvrnorm(n, mu, Sigma)
sim <- t(apply(x, 1, function(y) c(mean(y), sd(y))))
plot(sim, pch=16, cex=0.5, col="#00000008",
xlab="Mean", ylab="SD")
i <- order(sim[, 1])
lines(sim[i, 1], lowess(sim[i, 2], f=1/20)$y, col="Red", lwd=2)
# g <- cut(sim[, 1], quantile(sim[, 1], seq(0, 1, by=0.025)))
# boxplot(sim[, 2] ~ g)