/** # A Simple surface energy budget boundary condition For an atmospheric medium with a viscosity $\nu$, the buoyancy budget at the surface may read, $$G + B = Q_n,$$ With $G$, the soil flux, expressed using a lumped-parameter (i.e. a linear feedback) model: $$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, $$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, $$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](https://en.wikipedia.org/wiki/Robin_boundary_condition) for $b$ at the bottom surface, $$\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: $$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: $$\Pi_1 = \frac{B_0}{B_1} \approx -6,$$ $$\Pi_2 = TN \approx 2000 \rightarrow 500,$$ $$\Pi_4 = \frac{\sqrt{B_0T}}{\Lambda} \approx 5000 \rightarrow 2500,$$ $$Pr = \frac{\nu}{\kappa} \approx 1,$$ $$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)](https://doi.org/10.1175/JAS-D-19-0023.1) for details. ## Results ![The evolution of the buoyancy field and its horizontally-averaged profile](lumpedrob/output.mp4) ![The buoyancy field structure via $\mathrm{log}\left(\| \nabla b\| + 1\right)$](lumpedrob/db.mp4) */ #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; }