Solving two-dimensional heat balance equation

Question

I am trying to solve the heat balance equation for vanishing components of the velocity (Eqs. (1.1.17), (1.1.24), and (1.1.25) in here, neglecting temporal evolution):

$ (\partial_x^2+\partial_y^2)T(x,y) = a I(x,y) \\ (v_x(x,y) \partial_x + v_y(x,y) \partial_y) \omega(x,y) - b(\partial_x^2+\partial_y^2)\omega(x,y)=c \partial_xT(x,y) $

where $a,b,c$ are some (material) constants. The function

$I=\exp(-2(x^2+y^2)/0.02^2)$

I assume to be Gaussian with a certain width. I can solve the first equation using NDSolve, but the second gives me an error (Nonlinear coefficients are not supported in this version of NDSolve). This is my attempt:

Intens[x_, y_] := Exp[-2 (x^2 + y^2)/0.02^2]

K = 0.0257;
\[Alpha] = 1.2*10^-6;
\[Beta] = 3.67*10^-3;
g = 9.81;
\[Nu] = 1.568*10^-5;

xMin = -0.1;
xMax = 0.1;
yMin = -0.1;
yMax = 0.1;

pde = {
   D[T1[x, y], {x, 2}] + D[T1[x, y], {y, 2}] == \[Alpha]/
     K Intens[x, y],
   T1[xMin, y] == 0,
   T1[xMax, y] == 0,
   T1[x, yMin] == 0,
   T1[x, yMax] == 0
   };
soln = NDSolve[pde, T1, {x, xMin, xMax}, {y, yMin, yMax}];

(* The result looks correct *)
ContourPlot[
 Evaluate[T1[x, y] /. soln], {x, xMin, xMax}, {y, yMin, yMax}, 
 PlotLegends -> Automatic, PlotRange -> All]

(* Next step gives errors *)
T1hlp = soln[[1, 1, 2]];
pde = {
   vx[x, y]*D[w[x, y], x] + 
     vy[x, y]*
      D[w[x, y], 
       y] - \[Nu] (D[w[x, y], {x, 2}] + 
        D[w[x, y], {y, 2}]) == \[Beta] g D[T1hlp[x, y], x],
   D[\[CapitalPsi][x, y], {x, 2}] + 
     D[\[CapitalPsi][x, y], {y, 2}] == -w[x, y],
   vx[x, y] == D[\[CapitalPsi][x, y], y],
   vy[x, y] == -D[\[CapitalPsi][x, y], x]
   };
soln = NDSolve[
   pde, {w, \[CapitalPsi], vx, vy}, {x, xMin, xMax}, {y, yMin, yMax}];

Any advice?