Unexpected result from contour plot - where is the color gradient?

Question

When I plot a contour plot of my (in)famous code :P I get some unexpected results that I don't quite understand.

My code (the one that has NDSolve)that has been edited by several authors on mathematica.SE for ease of copying and running.

From the interpolating function polynomial that my NDSolve generates, I define a temperature T and density r and plot them as follows:

fac = 0.75;
T = 1 - (Bi*hSol[x, y, fac*TRup]/(Bi*hSol[x, y, fac*TRup] + 1));
a = 9.2*10^-4;
b = 4.5*10^-7;
r = 1.33/(1 + a *(T) + 
         b*(T)^2); (*Density as a function of temperature*)

Plot3D[
  1 - (Bi*hSol[x, y, fac*TRup]/(Bi*hSol[x, y, fac*TRup] + 1)), {x, 
  0, 
    L}, {y, 0, L},

  PlotRange -> {{0, L}, {0, L}, {0, 3.5}},
  BaseStyle -> {FontWeight -> "Bold", FontSize -> 18}
  ]



ContourPlot[
  Evaluate[\!\(
\*SubscriptBox[\(\[PartialD]\), \(x\)]r\) + \!\(
\*SubscriptBox[\(\[PartialD]\), \(\(y\)\(\ \)\)]r\), {x, 0, L}, {y, 
      0, L}],

  PlotRange -> {{0, L}, {0, L}, {0, 3.5}},
  BaseStyle -> {FontWeight -> "Plain", FontSize -> 18},
  ColorFunction -> Gray
  ]

The Plot3D plots the temperature and then a top view looks like this:

whilst the contour plot of density r, which is a function of T looks, bewilderingly, like this:

The little bumps as seen in the temperature plot aren't captured by my contour plot.

Now, I can't tell if I am missing any option in the contour plot. I assumed that I'd see a smooth gradient and not a binary black and white sudden transition in the contour plot of density.

Any suggestions? Thoughts? Recommendations?