How do I calculate $\sigma^2$ within the $\varepsilon = N(0, \sigma^2)$ notation The error term for an AR(1) model is assumed to be $\varepsilon = N(0, \sigma^2)$, a mean of 0. An example of a source that mentions this:

The AR(1) model can be written in intercept form,
$$z_t = \beta_0 + \beta_1 z_{t-1} + a_t$$
where $t = 1, \dots, T$, $\beta_0$ is the intercept coefficient, $\beta_1$ is the slope coefficient and $a_t$ is assumed to be $\text{NID}(0,\sigma_a^2)$.

Additionally, the error term is also defined as having a constant variance, in other words, the sigma squared isn't a variable but a constant.
My question is how can I calculate the sigma squared/variance of the error term.
I assume the variance isn't the same for each data set, that it's rather dependent on the data set and that I need to do some calculation with the data that I have to find the specific variance.

Your model has three parameters: $\beta_0$, $\beta_1$ and $\sigma^2$. Based on your data, you want to find values for these parameters that would in some sense bring your model closest to the data. This may sound like a minor point but we don't calculate parameters; we estimate them. To calculate would imply that 1) the result of this computation is the "true" value and 2) there is only one way to do it, neither of which is true. Actually, there are many different ways of estimating these parameters, and different estimators have different properties.
You can get estimates of the parameters of an AR(1) process by using ordinary least squares (OLS). Note that you effectively lose 1 observation in that case, since you don't have $z_0$. This is sometimes referred to as conditional sum of squares (CSS), i.e. conditional on the first observation. There are other ways to estimate these parameters where you don't lose an observation like this (e.g. you can use maximum likelihood estimation and the Kalman filter with exact diffuse initialization).
When using OLS, you first get estimates of $\beta_0$ and $\beta_1$, which you can then use to get the residuals:
$$\hat{\varepsilon}_t = z_t - \hat{\beta}_0 - \hat{\beta}_1 z_{t-1}$$
for $t = 2,...,T$. The usual estimate of $\sigma^2$ in a general OLS context where you have $N$ observations and $p$ regressors (including the constant if there is one) is:
$$\hat{\sigma}^2 = \frac{1}{N-p} \sum_{i=1}^N \hat{\varepsilon}_i^2$$
For your AR(1), that would translate to (we only have $T-1$ residuals, and two regressors including the constant):
$$\hat{\sigma}^2 = \frac{1}{T-3} \sum_{t=2}^T \hat{\varepsilon}_t^2$$
It is, essentially, the variance of the residuals $\hat{\varepsilon}_t$, which you might have expected would be a good estimate of the variance of the true errors $\varepsilon_t$. You may have expected this instead:
$$\tilde{\sigma}^2 = \frac{1}{T-1} \sum_{t=2}^T \hat{\varepsilon}_t^2$$
This is also an estimate of $\sigma^2$. The difference is that this estimate is actually too small on average: it tends to underestimate $\sigma^2$ because the residuals are computed from the same data from which we estimated the parameters, and are thus less variable than the "true" errors. The estimate that divides by $T-3$ is unbiased. On the other hand, the estimate that divides by $T-1$ is the maximum likelihood estimator and enjoys certain properties associated with that. In a large sample, the numerical difference between the two is small.
Whatever software you used to fit the regression likely also reports an estimate of $\sigma^2$ (or a related value named like "residual sum of squares", "residual standard error", "mean squared error", etc). It is likely to be the unbiased estimate, although it's possible that it isn't. If you wrote the software yourself, it would be fine to use the unbiased estimate. Otherwise, you should check documentation to find out what is used.
As a final note, in case that wasn't clear, $\sigma^2$ is the variance of the error term, not the variance of the process $z_t$. The variance of the process depends both on $\sigma^2$ and $\beta_1$. If the process is stationary, then its variance is $\frac{\sigma^2}{1-\beta_1^2}$.
