This problem was studied by Gerard Maze in Natural Density Distribution of Hermite Normal Forms of Integer Matrices (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}.$$$$\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 On random nonsingular Hermite Normal Form (2015, paywall). I quote their result:
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 $B$$M$, while Hu et al. choose the determinant uniformly from that interval.