It follows immediately from the special case (with $m=n$ and $s=t=k$) of Theorem 1 of Canfield and Mckay that the probability -- say $p_{n,k}$ -- in question is $$\sim\frac1{\sqrt e}\dfrac{\displaystyle{\binom{n}{k}^n}}{\displaystyle{\binom{n^2}{k n}}} \tag{1}\label{1}$$ if $n\to\infty$ and either $k=o(n^{1/2})$ or $k\asymp n-k\asymp n$. (Indeed, writing $M_{ij}:=1(j\in S_i)$ for all $i$ and $j$ in $[n]:=\{1,\dots,n\}$, we have a bijection between the set of all "even" $k$-coverings $(S_1,\dots,S_n)$ of the set $[n]$ and the set of all $0$-$1$ "incidence" matrices $M=(M_{ij}\colon(i,j)\in[n]^2)$ with all row sums and all column sums equal $k$.)
According to Conjecture 1 of Canfield and Mckay, \eqref{1} holds uniformly over all $k\in\{1,\dots,n-1\}$ as $n\to\infty$.
A prehistory of this problem is discussed in Canfield and Mckay as well.
For $k=2,3,4$, $$p_{n,k}=C_{n,k}\Big/\binom nk^n,$$ where $(C_{n,2})$, $(C_{n,3})$, and $(C_{n,3})$ are A001499, A001501, and A058528, respectively.