When we create vertical profiles (and hovmoller plots), we compute layer averages. Instead, we should define z-levels (probably the reference z-level profile from refZBot or whatever by default) and first interpolated to those levels before averaging horizontally (rather than over layers).
https://github.com/MPAS-Dev/MPAS-Analysis/blob/develop/mpas_analysis/ocean/ocean_regional_profiles.py#L456-L466