Pr = 0.71; Ra = 10^4;
Needs["NDSolve`FEM`"]
\[CapitalOmega] = ImplicitRegion[True, {{x, 0, 1}, {y, 0, 1}}];
op = {D[u[x, y], x]+D[v[x, y], y] \[Equal] 0, u[x, y]*D[u[x, y], x]+v[x, y]*D[u[x, y], y]+ D[P, x]-Pr*(D[u[x, y], x, x]+D[u[x, y], y, y]) \[Equal] 0, u[x, y]*D[v[x, y], x]+v[x, y]*D[v[x, y], y]+D[P, y] Pr*(D[v[x, y], x, x]+D[v[x, y], y, y])+Ra*Pr*\[Theta][x, y] \[Equal] 0, u[x, y]*D[\[Theta][x, y], x]+v[x, y]*D[\[Theta][x, y], y]-D[\[Theta][x, y], x, x]-D[\[Theta][x, y], y, y] \[Equal] 0}
bcs = {DirichletCondition[{u[x, y] \[Equal] 0, v[x, y] \[Equal] 0, \[Theta][x, y] \[Equal] 0}, x \[Equal] 0], DirichletCondition[{u[x, y] \[Equal] 0, v[x, y] \[Equal] 0, \[Theta][x, y] \[Equal] 1}, y \[Equal] 0], DirichletCondition[{u[x, y] \[Equal] 0, v[x, y] \[Equal] 0, \[Theta][x, y] \[Equal] 0}, x \[Equal] 1]}
NDSolveValue[{op, bcs}, {u, v, \[Theta], P}, {x, y} \[Element] \[CapitalOmega], Method -> {"FiniteElement"}]
