TLDR; $s(z)$ is asymptotically normal, and its variance is $\frac {12} {n-1}$ according to CLT for Markov chains. It can be shown that the distribution is a special case of generalized $\chi^2$ distribution.
Markov Chain approach, asymptotics and variance
The sequence $z_i$ is Markov chain because once you know $z_i$ the value of $z_{i+1}$ doesn't depend on $z_k$ where $k<i$. Therefore, Markov chain CLT is applicable. Here's how we apply it.
The sum, or any linear combination of normal r.v.s, is a r.v. itself. Knowing that $z_i\sim\mathcal N(0,2)$ or $z_i\sim\sqrt 2\space \mathcal N(0,1)$, we know that $z_i^2\sim 2\space\chi^2_1$, see the definition of $\chi^2$ distribution. Thus, $\sigma_z^2=\operatorname{var}[z_i^2]=2^2\times 2=8$.
Markov chain CLT states:
$$\sqrt{n-1}(s(z)-\mu)\sim\mathcal N(0,\sigma^2),$$
where $\mu=E[z_i^2]$ and $\sigma^2 = \sigma_z^2 + 2\sum_{k=1}^\infty \operatorname{cov}( z_{1}^2, z_{1+k}^2)=8+2\times 2=12$.
Hence $\operatorname{var}[s(z)]=\frac{12}{n-1}$
Here's the proof by simulation (Python):
import numpy as np
n = 51
s = np.mean(np.diff(np.random.randn(10000,n))**2,axis=1)
vars = np.var(s)
print(vars)
print(12/(n-1))
Output:
0.23526746023519335
0.24
Note, that if $z_i^2$ weren't correlated then $s(z)$ would have been from scaled $\chi^2$ distribution with variance $\frac 8 {n-1}$. However, due to overlapping terms $x_i$ in $z_i$ and $z_{i+1}$ we had to apply modified CLT to obtain the asymptotic distribution of $s(z)$.
Acknowledgements:
My initial answer, which I updated a few times already, did not account for correlation, which was pointed out by @Sextus Empiricus. Also, I used this answer for $\operatorname{cov}( z_{1}^2, z_{1+k}^2)$, where the correlation $\rho=\operatorname{corr}[z_i,z_{i+1}]=-1/2$ and we know that correlation disappears between $z_i$ and $z_j$ when $|i-j|>1$.
The distribution
Let's start with independent row vector of randoms $X'=(x_1,\dots,x_n)$. We get the row vector of differences $Z'$ aby pplying Toeplitz matrix $B'$ as follows $Z'=X'B'$, where
$$B' = \begin{bmatrix} {-1} & 0&\dots & 0 &0\\
1 & -1 & \dots&0 & 0 \\
0 & {1}& \dots & 0& 0\\
\vdots & \vdots & \vdots & \vdots & \vdots \\
0& 0 & \dots & 1 & -1 \\
0& 0 & \dots & 0 & {1}
\end{bmatrix}$$
Your quantity then is a quadratic form
$$s(z)=\frac 1 {n-1} X'B'BX$$
where $B'B$ has a form of a tridiagonal Toeplitz matrix:
Let's apply eigen decomposition $B'B=P'\Lambda P$ then we have:
$$s(z)=\frac 1 {n-1} X'P'\Lambda PX=\frac 1 {n-1} Y'\Lambda Y$$
where $Y=PX\sim\mathcal N(0,I_{n-1})$, i.e. each $Y_i$ (principal component) is an independent normal.
Hence, $$s(z)=\frac 1 {n-1} \sum_{i=1}^{n-1}\lambda_i Y_i^2$$
where $Y_i^2\sim\chi^2_1$ and $\lambda_i$ are eigenvalues. The eigenvalues of Toeplitz tridiagonal matrices are known to form a sine wave and easy to find, see "The Eigenproblem of a Tridiagonal P-Toeplitz Matrix" by Gover .
So the distribution can be seen as a linear combination of $\chi^2$ variables or a generalized $\chi^2$ distribution.
Miscellaneous
We can define a row vector $V'=(x_1,z_1,\dots,z_{n-1})$, then it can be obtained by applying a matrix $D'$ to the original observations $V'=X'D'$, the matrix $B'$ above is the subset of columns of $D'$:
The matrix $D'D$ looks like this:
We can get the matrix $U'$ that recovers the original vector from $V$ as follows: $X'=V'U'$, and $U'=D'^{-1}$.
The matrix $U'$ is upper unit triangular, meaning $u_{ij}=1_{i\ge j}$:
Matrix $A=U'U$, which appears in the quadratic form, has a very interesting form: $a_{ij} = n+1-min(i,j)$, e.g. $n=5$:
