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}.$$

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). I quote their result (since the paper is behind a paywall):

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