lumpedrob.c

A Simple surface energy budget boundary condition
- Results

A Simple surface energy budget boundary condition

For an atmospheric medium with a viscosity \nu, the buoyancy budget at the surface may read,

\displaystyle G + B = Q_n,

With G, the soil flux, expressed using a lumped-parameter (i.e. a linear feedback) model:

\displaystyle G = \Lambda (b_{\mathrm{surf}} - b_0).

with b_0 a reference buoyancy value and \Lambda the lumped soil conduction parameter. For the sensible buoyancy flux B, per definition,

\displaystyle B = \kappa \frac{\partial b}{\partial \mathbf{n}},

with \kappa the heat conductivity of the atmospheric medium. Finally, for a simple diurnal cycle with period T, Q_n is written as,

\displaystyle Q_n = \mathrm{max} \left[ B_0\ \mathrm{sin}\left( \frac{2\pi t}{T} \right),\ B_1 \right],

with B_0 and B_1 buoyancy flux scales. Giving a Robin boundary condtion for b at the bottom surface,

\displaystyle \Lambda b_{\mathrm{surf}} + \kappa \frac{\partial b}{\partial \mathbf{n}} = Q_n + \Lambda b_0.

Initially, the atmospheric fluid is at rest with a constant (vertical y) stable stratification:

\displaystyle b(t=0) = b_0 + N^2y,

From the 7 system parameters (in order of appearance: \{ \nu, \Lambda, \kappa, B_0, T, B_1, N\}), we can construct 5 dimensionless parameters:

\displaystyle \Pi_1 = \frac{B_0}{B_1} \approx -6, \displaystyle \Pi_2 = TN \approx 2000 \rightarrow 500, \displaystyle \Pi_4 = \frac{\sqrt{B_0T}}{\Lambda} \approx 5000 \rightarrow 2500, \displaystyle Pr = \frac{\nu}{\kappa} \approx 1, \displaystyle Re = \frac{B_0}{\nu} \left( \frac{2T}{\pi N^2} \right) ^{2/3} \approx 10^8 \rightarrow 5000.

See Van Hooft et al (2019) for details.

Results

The evolution of the buoyancy field and its horizontally-averaged profile

The buoyancy field structure via \mathrm{log}\left(\| \nabla b\| + 1\right)

#include "navier-stokes/centered.h"
#include "tracer.h"
#include "diffusion.h"

#define robin(a,b,c) ((dirichlet ((c)*Delta/(2*(b) + (a)*Delta)))	\
		      + ((neumann (0))*					\
			 ((2*(b) - (a)*Delta)/(2*(b) + (a)*Delta) + 1.)))

#define QN (max(B0*sin(2*pi*t/T), B1))

double T = 24*60, Lc = 1;
double Pi1 = -6, Pi2 = 500, Pi4 = 2500, Re = 5000;
double LAMBDA, B0, B1, sqN, NU;

int maxlevel = 8;

scalar b[], * tracers = {b};
b[bottom] = robin (LAMBDA, NU, QN);
b[top] = dirichlet (sqN*y);

FILE * gp;
face vector av[];

int main() {
  sqN    = sq(Pi2/T);
  B0 = sq(Lc)*sqN*pi/(2.*T);
  B1     = B0/Pi1;
  LAMBDA = sqrt(B0*T)/Pi4;
  NU = B0*pow(2*T/(pi*sqN), 2./3.)/Re;
  L0     = 3*Lc;
  X0     = -L0/2;
  periodic (left);
  u.t[bottom] = dirichlet (0.);
#if (dimension == 3)
  u.r[bottom] = dirichlet (0.);
  periodic (back);
#endif
  const face vector muc[] = {NU, NU, NU};
  mu = muc;
  a = av;
  run();
}

event init (t = 0) {
  DT = 1./(sqrt(sqN)*10.);
  foreach() {
    b[] = sqN*y;
    foreach_dimension()
      u.x[] = noise()/1000.;
  }
  boundary ({b, u});
  gp = popen ("gnuplot", "w");
  fprintf(gp,
	  "set term pngcairo size 500, 500\n"
	  "set xr [%g: %g]\n"
	  "set yr [0 : %g]\n"
	  "set grid\n"
	  "set size square"
	  "set key topleft"
	  "set xlabel 'b'\n"
	  "set ylabel 'y/L_c'\n",
	  B1/LAMBDA, sqN*Lc*2., Lc*2.);
}

event acceleration (i++) {
  coord grav = {0, 1, 0};
  boundary({b});
  foreach_face()
    av.x[] = grav.x*face_value(b, 0);
}

event tracer_diffusion(i++) 
  diffusion (b, dt, mu);

event adapt (i++) {
  double be = sqN*Lc/40., ue = 0.05;
  adapt_wavelet ({b, u}, (double[]){be, ue, ue, ue}, maxlevel, 4);
}

event damp (i++) {
  foreach() {
    if (y > 2.*L0/3.) {
      foreach_dimension()
	u.x[] *= exp (-dt*(y - 2*L0/3));
      b[] -= (b[] - sqN*y)*(1. - exp (-dt*(y - 2*L0/3)));
    }
  }
  boundary ({b, u});
}

#include "profile5c.h"
int frame = 0;
event movies (t += T/1000.) {
  scalar db[];
  foreach() {
    db[] = 0;
    foreach_dimension() 
      db[] += sq((b[1] - b[-1])/(2*Delta));
    db[] = db[] > 0 ? log(sqrt(db[]) + 1.) : 0;
  }
  boundary ({db});
  output_ppm (b, file = "b.mp4", linear = true,
	      n = 500, min = 0, max = sqN*Lc);
  output_ppm (db, file = "db.mp4", linear = true,
	      n = 500, min = 0, max = 3*log(sqN + 1));
  double by, yp = 1e-6;
  fprintf (gp,
	   "set output 'plot%d.png'\n"
	   "set title 'Time since sunrise:  %02d:%02d (HH:MM)'\n"	
	   "plot '-' w l lw 3 t 'buoyancy'\n",
	   frame, (int)(t/60.), (int)floor(fmod(t, 60.)));
  while (yp < L0) {
    by = 0;
    double dy = average_over_yp ({b}, &by, yp);
    fprintf (gp,"%g %g\n", by, yp);
    yp += dy;
  }
  fprintf (gp,"e\n");
  frame++;
}

event stop (t = 2.*T) {
  pclose (gp);
  system ("rm mov.mp4");
  system ("ffmpeg -f image2 -i plot%d.png -c:v libx264 -vf format=yuv420p -y mov.mp4");
  system ("rm plot*");
  sleep (10);
  system ("ffmpeg -y -i b.mp4 -i mov.mp4 -filter_complex hstack output.mp4");
  return 1;
}