This problem was studied by Gerard Maze in <A HREF="https://arxiv.org/abs/1009.4826">Natural Density Distribution of Hermite Normal Forms of Integer Matrices</A> (2010). One result:

The probability that an $n\times n$ integer matrix $A$ has an Hermite normal form (HNF) with the integers $d_1,d_2,\ldots,d_{n−1},d$ on the diagonal is given by
$$\bigl(\zeta(n)\zeta(n-1)\cdots\zeta(2)d_1^n d_2^{n-1}\cdots d_{n-1}^2\bigr)^{-1},$$
with $\zeta(n)$ the Riemann zeta function.

The expectation value and variance of a diagonal element have been calculated by Hu et al. in <A HREF="https://www.sciencedirect.com/science/article/pii/S0022314X16000354">On random nonsingular Hermite Normal Form</A> (2015, paywall). I quote their result:

<IMG SRC="https://ilorentz.org/beenakker/MO/HNF.png"/>

Note that the distribution of the matrices is different in the two papers: Maze chooses the matrix elements uniformly in the interval $(-M,M)$ for large $M$, while Hu et al. choose the determinant uniformly from that interval.