Q - vector in the horizontal plane

Definition of Q1 and Q2

Q1[x_,y_,Θ_,Ψ_]:= (∂_(x,y))Ψ[x,y](∂_x)Θ[x,y]-(∂_(x,x))Ψ[x,y](∂_y)Θ[x,y];
Q2[x_,y_,Θ_,Ψ_]:= (∂_(y,y))Ψ[x,y](∂_x)Θ[x,y]-(∂_(x,y))Ψ[x,y](∂_y)Θ[x,y];

The geostrophic velocity is here described byits streamfunction Ψ. The temperature is noted as θ. 
It is known that θ=(∂_z)Ψ 
up to some scaling factors depending on the choices of variables and vertical coordinates. However, we treat Ψ and θ as two independant functions in the horizontal (isobaric) plane.
The scales are here arbitrary.

Plot module

First graph
Temperature : cold (blue) to warm (red)
Geostrophic streamfunction or geopotential : black lines : positive values (solid), negative values (dashed)
(OverVector[Q]) - vector : vectors
Second graph
Geostrophic streamfunction and divergence of the (OverVector[Q]) - vector:  positive values (blue), negative values (purple)
Convergence of (OverVector[Q]) is associated with ascending motion while divergence of (OverVector[Q]) is associated with convergence. Attention: it is necessary to known the 3D distribution of (OverVector[Q]) and to solve the Poisson equation -2(OverVector[∇]). (OverVector[Q]).
The scales on the axis of the graphs are arbitrary. The mean flow, when it is relevant, is always oriented left to right.

CoolColor[ z_ ] :=RGBColor[1-z,z,1];
Qdisp[θ_,Ψ_,segx_,segy_]:=Module[{g1,g2,g3},
g1=VectorPlot[Evaluate[{Q1[x,y,θ,Ψ],Q2[x,y,θ,Ψ]}],Evaluate[segx],Evaluate[segy],AspectRatio->0.5,VectorStyle->GrayLevel[0.4]];
g2=ContourPlot[θ[x,y],Evaluate[segx],Evaluate[segy],ContourStyle->None,Contours->50,ColorFunction->(RGBColor[Min[2 ((#1)^2),1],4 #1 (1-#1),Min[2 (((1-#1))^2),1]]&),AspectRatio->0.5];
g4=ContourPlot[UnitStep[Ψ[x,y]]Ψ[x,y],Evaluate[segx],Evaluate[segy],ContourShading->None,ContourStyle->Black,Contours->5,AspectRatio->0.5];
g5=ContourPlot[(1-UnitStep[Ψ[x,y]])Ψ[x,y],Evaluate[segx],Evaluate[segy],ContourShading->None,ContourStyle->{{Black,Dashed}},Contours->5,AspectRatio->0.5];
g3=ContourPlot[Evaluate[(∂_x)Q1[x,y,θ,Ψ]+(∂_y)Q2[x,y,θ,Ψ]],Evaluate[segx],Evaluate[segy],AspectRatio->0.5,ColorFunction->CoolColor];
GraphicsColumn[{Show[g2,g4,g5,g1],Show[g3,g4,g5]}]]

Case of the surface wave

θa=Function[{x,y},-y];
Ψa=Function[{x,y},((Cos[x]Cos[y])/(Cosh[1.5 y]))];
segxa={x,-π,2 π};segya={y,-π/2,π/2};

Qdisp[θa,Ψa,segxa,segya]

Case of confluence

θb=Function[{x,y},-y];
Ψb=Function[{x,y},-x Tanh[1.5 y]- y];
segxb={x,0,3π};segyb={y,-π/2,π/2};

Qdisp[θb,Ψb,segxb,segyb]

Case of jet streak

θc=Function[{x,y},-y];
Ψc=Function[{x,y},(((x Tanh[x] -2π))/(Cosh[y]))Tanh[1.5 y]- 2 y];
segxc={x,-1.5π,1.5π};segyc={y,-π/2,π/2};

Qdisp[θc,Ψc,segxc,segyc]

Case of jet separation

θd=Function[{x,y},-y];
Ψd=Function[{x,y},-x Tanh[x] ((Tanh[1.5 y])/(Cosh[y]))- 2 y];
segxd={x,-1.5π,1.5π};segyd={y,-π/2,π/2};

Qdisp[θd,Ψd,segxd,segyd]

Warm advection

θe=Function[{x,y},-y-0.5 x];
Ψe=Function[{x,y},-Tanh[1.5 y]- 2 y];
segxe={x,-1.5π,1.5π};segye={y,-π/2,π/2};

Qdisp[θe,Ψe,segxe,segye]

Cold advection

θf=Function[{x,y},-y+0.5 x];
Ψf=Function[{x,y},-Tanh[1.5 y]- 2 y];
segxf={x,-1.5π,1.5π};segyf={y,-π/2,π/2};

Qdisp[θf,Ψf,segxf,segyf]

Waves of temperature and geopotential

No variation in latitude is assumed.
Notice that the (OverVector[Q]) vector does not depend of the phase of temperature with respect to the geopotential.

Geopotential in phase with temperature

θg=Function[{x,y},-2y+Cos[x]];
Ψg=Function[{x,y},-  y+ Cos[x]];
segxg={x,-1.5π,1.5π};segyg={y,-π/2,π/2};

Qdisp[θg,Ψg,segxg,segyg]

Geopotential in quadrature with temperature

θh=Function[{x,y},-0.5 y+Sin[x]];
Ψh=Function[{x,y},- y+ Cos[x]];
segxh={x,-1.5π,1.5π};segyh={y,-π/2,π/2};

Qdisp[θh,Ψh,segxh,segyh]

Geopotential in opposition with temperature

θi=Function[{x,y},-2y-Cos[x]];
Ψi=Function[{x,y},-  y+ Cos[x]];
segxi={x,-1.5π,1.5π};segyi={y,-π/2,π/2};

Qdisp[θi,Ψi,segxi,segyi]

Sandbox

θz=Function[{x,y},-y];
Ψz=Function[{x,y},-y - (y^2)+((Cos[x-2 y]Cos[y])/(Cosh[1.5 y]))];
segxz={x,-π,2 π};segyz={y,-π/2,π/2};

Qdisp[θz,Ψz,segxz,segyz]