Pressure on the surface for case with free surface boundary condition

Hello folks,

My case involves free surface boundary condition, so I turned off Pressure normalization. I am curious about that the pressure are not negative everywhere on the surface when I prescribed more than 1 crustal domains with different crust thickness and effective viscosity (but with the same density; lower viscosity mantle is below the crust). The region with thinner but higher viscosity crust has positive pressure and subsidence on the surface. Uplift is found on the surface of other regions. However, if only one crust domain is considered, the pressure on the surface is negative all the time, which is what I expect.

I guess the issue lies in the mesh refinement, for example, inaccuracies of meshing the domain that was introduced by an unreasonable set to “set Strategy”. However, it’s not the reason.

Does anybody run into a similar issue?

Best,
Shaozhuo