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}.$$
The expectation value and variance of a diagonal element have been calculated by Hu et al. in On random nonsingular Hermite Normal Form (2015). I quote their result (since the paper is behind a paywall):
