Is there a way to have heatmap and contours from other variables in Plots.jl with GR?

Specific Domains Visualization
,
2

You could shift the additional datasets and their contour levels to the mean z0 of the primary dataset and request that the z0 contour be plotted for each shifted version:

z0 = mean(mag)

p = heatmap(P_range, δ_range, mag; xlabel="P", ylabel="δ", title="‖∇Q_T‖", color=:viridis, clims=cl)

levels1 = [-0.15, 0.15]
for lvli in levels1
    contour!(P_range, δ_range, z0 .+ detJ .- lvli; levels=[z0], clims=cl, lc=:black, ls=:dash, lw=1.5)
end

levels2 = [0.0]
for lvli in levels2
    contour!(P_range, δ_range, z0 .+ Q_total .- lvli; clims=cl, levels=[z0], lc=:green, lw=2)
end

p