That is interesting. Thanks for pointing it out. In the range of operation, e.g. 0 <= z <= 1 and 140 < p < 180 (pressure limits are arbitrary), the function dzRHS is continuous. Demonstrated by Plot3D[dzRHS[z, p], {z, 0, 1}, {p, 140, 180}].

It could be that z is getting outside of this region, in this example z might be getting slightly above one. The intent of the piecewise implementation of dzRHS was to keep z between 0 and 1. Maybe this approach isn't successful at limiting z.

Is it possible that with the WhenEvent implementation, NDSolve looks more carefully at the boundary (z near 1 in this example).

Your thoughts?