Just in case, I'll outline a quantitative approach to the problem. Just for the record, as has been already mentioned many times, fig. 2 is the correct picture.

I will consider a uniformly, axially magnetised ring, which should be an accurate model for say a neodymium magnet. It is the same as two concentric solenoids of opposing currents with uniform surface currents. It is also equivalent to two stacked annuli that are oppositely and uniformly magnetically charged. This gives you two approaches to compute the magnetic field. Technically, the charge approach gives you $H$, but outside the magnet, this is the same as $B$. The former is physically motivated while the latter is conceptually more convenient.

To make things quantitative I'll use cylindrical coordinates about the ring's axis. Say your ring has inner radius $R_i$ and outer radius $R_o$ and height $D$ and magnetisation $\vec M = M_z\vec e_z$, then the outer and inner surface current density is (using cylindrical coordinates):
$$
\begin{align}
\vec J_o &= M_z \vec e_\phi & \vec J_i &= -M_z \vec e_\phi
\end{align}
$$
while the top and bottom charge density is:
$$
\begin{align}
\sigma_t &= -M_z & \sigma_b &= M_z
\end{align}
$$
In your examples $M_z>0$. The outer surface current rotates counterclockwise and the inner surface rotates clockwise (when viewed from above). The top face is therefore positively charged and the bottom face negatively charged.

This gives you two approaches to build intuition on the qualitative behaviour of the field lines and to quantititatively calculate them. The exact formula has no simple expression and involves double integrals. It gives the same result regardless whether you use Biot-Savart's law following the current approach or Coulomb's law following the charge approach.

Both approaches explain why there is a "polarity inversion" in the "air gap" (the hole of the ring). From the current perspective, the outside surface currents generate an upward field in the gap, but the inside surface currents generate a downward field (the usual solenoid field). It is the latter that dominates after superposition since it is closer. Indeed the field decays with distance to the source, and this can be made rigorous on the axis. You can reach the same conclusion for the entire air gap with the charge perspective. You have positive charge above the gap and negative charge below, so the field is necessarily directed downwards. In either case, the result is mathematically airtight (without having to calculate exactly the field): the field points downwards in the gap.

This is consistent with your equation of $B$ on the $z$ axis. No, the formula is not meant to only be applied above the top face, it is valid on the entire axis. In fact, using the charge approach, you can trace back where each term comes from. Each annulus can be viewed as overlapping oppositely, uniformly charged disks of radius $R_i,R_o$. The four terms come from each of the four disks. Just as Coulomb's law is valid everywhere in space, this formula is valid everywhere on the axis.

From the current/magnetisation perspective, it's a bit trickier to interpret the terms. Naturally the first two terms correspond to the outer solenoid and the next bracketed to terms to the inner solenoid. In terms of magnetisation, you need to imagine the ring to be the superposition of two concentric oppositely, uniformly magnetised cylinders of respective radius $R_o,R_i$. This is why they correspond to the previous formula in your linked page. You cannot understand where the two sub terms come from as from the current/magnetization perspective, you cannot have a net magnetic charge. You can almost get there by adding and subtracting by the same constant. Essentially, a finite solenoid can be viewed as two semi infinite solenoids with an overlap and an infinite solenoid to cancel the non overlapping region. Therefore, up to an additive constant, each of the four terms can be viewed as the field of an infinite axially magnetised cylinder with a discontinuity at a certain cross section, or equivalently, an infinite solenoid with a discontinuity of the surface current at a cross section.

The general result can be plotted with any existing visualisation tool. They all confirm your "polarity flip." This backed by both of the ressources you shared: [falstad](https://www.falstad.com/vector3dm-old/) and [R&J magnetics](https://tinyurl.com/bdduxmpb).

Hope this helps.