The units of the NDF are tricky. For whatever it's worth, Heitz's convention of defining it relative to a 1 m² reference geometric surface is unusual, and although I can see why he would want to define it that way for conceptual simplicity, it does not really match how NDFs are used in practice. I definitely had a head-tilt moment when I first read that passage in his paper.
The problem with NDFs is that they involve two different area measures: a macrosurface (or "geometric") area and a microsurface area. As far as I understand, the true and complete definition of the units of the NDF is:
$$
\frac{\text{micro area}}{\text{solid angle} \cdot \text{macro area}}
$$
In other words, if we double-integrate the NDF over a region of solid angle $\Omega$ and a region of macrosurface $A$, we obtain the area of the microsurface that has normals within that region of solid angle:
$$
\int_\Omega \int_A D(\omega)\,dA\,d\omega = \text{total area of microfacets above $A$ with normals in $\Omega$}
$$
Note that this is dimensionally consistent, even maintaining the distinction between the two area measures. The $dA$ is an element of macro area, and $d\omega$ an element of solid angle, so those cancel out the denominator of $D(\omega)$. The numerator is then micro area, which matches the right hand side.
Now if we forget the distinction between these two area measures and cancel the area units against each other, we obtain just $1/\text{sr}$, which is consistent with the vast majority of literature. But there is a subtlety here. Usually, when we see a distribution that has the units of $1/\text{sr}$, this means it's a distribution over the sphere and that if you integrate it over the sphere you should get 1. But note that if you integrate the NDF over the sphere, you don't in general get 1. Instead, as the above definition shows, you get the overall ratio of micro area to macro area. For a rough surface, this is greater than 1 (and the rougher it is, the larger the value), as the microsurface is crinkly and rough and therefore has more area than the corresponding macrosurface.
The correct normalization condition for an NDF is
$$
\int_{S^2} D(\omega) (n \cdot \omega) \, d\omega = 1
$$
(as can be readily verified if you take any of the usual NDF formulas, set a nonzero roughness, and try numerically integrating it—they don't come out to 1 unless you use the above integrand).
The reason for the $(n \cdot \omega)$ factor in the integrand is to convert microsurface area into the corresponding (projected) macrosurface area. In other words, this factor has the units:
$$
n \cdot \omega = \frac{\text{macro area}}{\text{micro area}}
$$
With this, the area units of the NDF really do cancel properly and you're left with a true spherical distribution that integrates to 1.
Now why does Heitz say something different? He wants to eliminate the "macro area" factor by setting it to a standardized value of 1 m², then defining the NDF as simply $\text{micro area}/\text{solid angle}$. Thus its units become $\text{m}^2 / \text{sr}$, and the normalization condition becomes
$$
\int_{S^2} D(\omega) (n \cdot \omega) \, d\omega = 1\ \text{m}^2
$$
(cf. equation (9) in the paper). It's no longer a dimensionless scalar 1, but an area of 1 square meter.
I don't love this definition. On the one hand, it does make it really explicit that an NDF involves area "somehow", and isn't simply a spherical distribution. It avoids the somewhat bizarre situation of having two area measures whose units cancel, but whose values don't cancel. On the other hand, in real life we obviously don't apply NDFs only to flat surfaces of 1 m². We want to apply NDFs to surface patches of any size we like—such as the patch enclosed in a pixel when rendering. The obvious way to do that is by normalizing out the 1 m², so that the NDF is expressed as a fractional area relative to any desired patch size. And then you get back to the definition I stated at the beginning of this answer.
Indeed, even Heitz isn't really consistent with this choice, as you noted. Starting from equation (14) in the paper, he drops the implicit 1 m² factors that ought to be there. When he proceeds to define the distribution of visible normals, $D_{\omega_o}(\omega_m)$, he no longer says that the units are $\text{m}^2/\text{sr}$, but just the usual $1/\text{sr}$ (cf. Table 2 in the paper).
To be fair, Heitz's distribution of visible normals is defined in such a way that it is a legitimate spherical distribution by itself, with the usual normalization condition:
$$
\int_{S^2} D_{\omega_o}(\omega_m) \, d\omega_m = 1
$$
(cf. equation (18) in the paper), so it is no longer necessary to think about some hidden area units in order to understand the meaning of this function.
Arguably BRDFs and NDFs should have been defined this way as well all along, incorporating the $(n \cdot \omega)$ factor into the function, so that they would be legitimate spherical distributions and we could avoid all this hassle. Sadly, due to historical path-dependence we didn't end up there.
So, to sum up, the oddities of the NDF have their origin in the (IMO under-appreciated and under-explained) fact that it has two different area measures hiding in it. If you just look at the units, they appear to cancel out, but nevertheless you need to keep the areas in mind in order to understand NDFs properly.
