canopy.c

Tempeature-gradient sharpening in a canopy.

Tempeature-gradient sharpening in a canopy.

The evolution of a temperature profile is studied with the minimal(?) ingredients for stable layer formation. Consider a cool bottom surface below a warmer canopy top (z = z_c) with a fixed difference \Delta T. The evolution of the temperature profile is described with an effective diffusivity K,

\displaystyle \frac{\partial T}{\partial t} = \frac{\partial}{\partial z}K\frac{\partial T}{\partial z }.

K is modelled by considering 3 effects:

A constant Background diffusivity due to shear, top down mixing etc (K_b).
A Height-dependend diffusivity correction for the canopy details.
A stability correction.

In a dimenonless setting, we choose z_c = 1, \Delta T = 1 and K_b = 1. The canopy correction defines K_c:

\displaystyle K_c = K_b \left( 1 + 4\frac{C_k}{z_c^2} \left(z\left(z - z_c\right)\right)\right),

#define Kc (Kb*(1 + 4*Ck*(x*(x - 1))))

Which gives a maximum of K_b at top and bottom and a minimum of (1-C_k)K_b halfway the canopy. This is to model the effects of the turbulent air induced by the surface of the Earth and the canopy crown.

The stability correction gives K_s,

\displaystyle K_s = K_c e^{-a\frac{\partial T}{\partial z}},

with a an inverse temperature-gradient scale.

Finally, to omit the most prominent issue with this description, there is an minimum value for K; K_{\mathrm{min}}.

\displaystyle K = \mathrm{max} \left(K_s, K_{\mathrm{min}}\right).

#define K(grad) (max(Kc*exp(-a*(grad)), Kmin))

Numerical set-up

The system is set up and solved for

#include "grid/multigrid1D.h" // a 1D grid
#include "diffusion.h"         
#include "run.h"              // Timeloop

The values of the physical parameters are chosen adhoc to produce a layered scenario.

double Kb = 1, Kmin = 0.025, Ck = 0.4, a = 0.6;

double tend = 15;

scalar T[];
T[left]  = dirichlet (0);  //bottom
T[right] = dirichlet (1);  //top

int main() {
  N  = 512;
  DT = 0.01;
  run();
}

event init (t = 0) {
  foreach()
    T[] = x;                // Linearly stratified
  boundary ({T});
}

event diff (i++) {
  dt = dtnext (DT);
  face vector D[];
  foreach_face() 
    D.x[] = K (face_gradient_x (T, 0));
  boundary ({D.x});
  diffusion (T, dt, D);
}

Output

A layer emerges

The movie is generated with gnuplot and ffmpeg.

FILE * gnuplotPipe;

event init (t = 0) {
  gnuplotPipe = popen ("gnuplot", "w");
  fprintf(gnuplotPipe,
	  "set term pngcairo\n"
	  "set xr [0: 1]\n"
	  "set yr [0: 1]\n"
	  "set key top left\n"
	  "set grid\n"
	  "set title 'Temperature and its Diffusivity'\n"
	  "set xlabel 'T/ΔT, K/K_b'\n"
	  "set ylabel 'z/z_c'\n");
}

int frame = 0;
event movie (t += 0.1) {
  fprintf (gnuplotPipe,
	   "set output 'plot%d.png'\n"
	   "plot 0.6*(x-0.5) + 0.5 w l lw 1 t 'Max. gradient', ",
	   frame);
  fprintf (gnuplotPipe, "'-' w l lw 5 t 'T', '' w l lw 4 t 'K'\n");
  foreach()
    fprintf (gnuplotPipe, "%g %g\n", T[], x);
  fprintf (gnuplotPipe, "e\n");
  foreach_face()
    fprintf (gnuplotPipe, "%g %g\n", K (face_gradient_x(T, 0)), x);
  fprintf (gnuplotPipe, "e\n");
  frame++;
}

event output_bart (t += 1) {
  char fname[99];
  sprintf (fname, "data%d", (int)(t + 0.5));
  FILE * fp = fopen(fname, "w");
  fputs("#Height\tT\tK\n", fp);
  foreach_face()
    fprintf (fp, "%g\t%g\t%g\n",x, face_value(T, 0),
	     K (face_gradient_x (T, 0)));
  fclose (fp);
}

event stop (t = tend) {
  system ("rm mov.mp4");
  system ("cp plot145.png pltend.png");
  system ("ffmpeg -r 25 -f image2 -i plot%d.png -c:v libx264 \
          -vf format=yuv420p -y mov.mp4");
  system ("rm plot*");
  return 1;
}

A criterion for layer formation

A steady state solution is characterized by a constant flux:

\displaystyle K\frac{\partial T}{\partial z} = H

For convinience we write,

\displaystyle y = \frac{\partial T}{\partial z}

then,

\displaystyle K_c e^{-\alpha y}y = H,

There exists a maximum for H when y = \alpha^{-1}. At this value for y the flux decreases with increasing gradients. It is well known that such systems Collapse (van de Wiel et al., van de Wiel et al.). Thus we have a constraint y < \frac{1}{\alpha}.

\displaystyle H_{\mathrm{max}} = \frac{K_c}{e\alpha},

\displaystyle H_{\mathrm{max}} = \frac{K_b \left(1-\frac{4}{z_c^2}C_k(z^2 - z)\right)}{e\alpha}.

The most sustainable heatflux is the minimal value of H_{\mathrm{max}}(z), which is at height z = 0.5z_c,

\displaystyle \frac{H_{\mathrm{max, s}}}{K_b} = \frac{1 - C_k}{e\alpha}.

Which would conclude the analysis if the system had flux boundary conditions. In this setup however, the steady-state heatflux H is a function of the parameters. To find the relation between the critical vales for C_k and \alpha, we need to invert for the profile of y.

\displaystyle y(z) = \frac{W\left(\frac{\frac{H}{K_b}}{1 - \frac{4}{z_c}C_k(z^2 - z)}\right)}{-\alpha}

where W(x) is the Lambert W function. Then we can then find H such that \langle y \rangle = \int_0^{z_c} y \mathrm{d}z = \frac{\Delta T}{z_c}. This will be a tedious excersize.

Side note: because the height-average of y is equal to \langle y \rangle = \frac{\Delta T}{z_c}, and y is maximum at z = 0.5z_c, we can write,

\displaystyle y_{z = 0.5z_c} > \frac{\Delta T}{z_c},

So if \alpha > \frac{z_c}{\Delta T}, the system is always unstable and collapes for any C_k.