Hello,

The desired system is an Mercator projection.

The Mercator projection is a map projection, without a defined height component. This allows you to divide your problem into two independent adjustment tasks. Separate your observations into planar observations such as GNSS-2D baselines, which should be projected into the target datum, direction and horizontal distances and adjust the horizontal network. For the height component, adjust a 1D network using the leveling data (and the height differences derived from the zenith angle and the slope distance, and - if required - the z-components of the baseline).

Seeing as the height is not referenced to the ellipsoid but rather local orthometric.

That's right. However, the ellipsoidal height is approximately the sum of the orthometric height and the geoid height, i.e. h = H + N. The difference is

Δh = hA - hb = HA + NA - HB - NB ≈ HA - HB,

if NA ≈ NB is assumed, which may valid in small networks. If NA differs significantly from NB, you must correct the geoid heights to combine these heights with your orthometric heights.

Kind regards
Micha