sandbox/Antoonvh/subsiding-shell.c

    At a cloud interface evaporation occurs due to mixing of moisture and temperature. Image via wideopen spaces.

    At a cloud interface evaporation occurs due to mixing of moisture and temperature. Image via wideopen spaces.

    Subsiding Shells at the edge of a cumulus-cloud’s updraft

    It is well known that mixing at a cumulus cloud’s interface can cause a so-called subsiding shell. Nair et al. (2019) propose a scenario to study these physics in an idealized setting. We follow their setup for the flow configuration.

    #include "navier-stokes/centered.h"
    #include "tracer.h"
    #include "diffusion.h"
    #include "fractions.h" //For the heavyside function;

    A sequence of DEFINEs is used to compute the cloud-liquid water content as a function of pressure and temperature. Note that there is no field dedicated to store the liquid water specific humidity, its prognostic.

    double Rv = 461.5;        //in J/kg/K
    double Rd = 287.;         //in J/kg/K
    double L = 2500000.;      //in J/kg
    double cpd = 1004.;       //in J/Kg/Kq
    #define P00 100000.       //A constant reference pressure
    double P0 = 101780., Pe = 7500; //Surface pressure and 1/e length scale
    
    #ifndef RHO0            
    #define RHO0 1.3 //Base density in kg/m^3  
    #endif
    #define pres  (P0*(exp(-1500./Pe)))  
    #define EXNER (pow(pres/P00, Rd/cpd))
    #define TK1   (thl[]*EXNER)
    #define ES    (610.78*exp(17.27*(TK1 - 273.16)/(TK1 - 35.86)))
    #define QSL   (Rd/Rv*(ES/(pres - (1. - Rd/Rv)*ES)))
    #define QSAT  (QSL*((1. + ((sq(L)/(Rv*cpd*sq(TK1)))*qt[]))/(1 + ((sq(L)/(Rv*cpd*sq(TK1)))*QSL))))
    #define QC    (max(qt[] - QSAT, 0))

    The virtual potential temperature in a cell (\theta_v,THV) can be computed from \theta_l, q_c and q_t.

    #define THETA (thl[] + L*QC/(cpd*EXNER)) 
    #define THV (THETA*(1. - (1. - Rv/Rd)*(qt[] - qt0) - Rv/Rd*QC))

    The virtual potential temperature is related to the buoyancy,

    double g_acc = 9.81;
    #define B_CENTERED (g_acc*((THV - T_ref)/T_ref))

    The scenario

    The case is defined by in-cloud and environment variables:

    double T_ref = 294., qt0 = 0.004;
    double qtc = 0.009, thvc = 289;
    double Lc = 2, wc = 1;
    scalar qt[], thl[], * tracers = {qt, thl};
    scalar H[]; //Heavyside function
    
    face vector av[];
    double tau; 
    int main() {
      H.prolongation = H.refine = fraction_refine; //Heavyside fraction
      periodic (bottom);
      L0 = 30;
      a = av;
      tau = Lc/wc;
      const face vector muc[] = {1./1000., 1./1000., 1./1000.};
      mu = muc;
      N = 128;
      run();
    }
    
    event init (t = 0) {
      fraction (H, Lc - x); //heavyside fraction 
      foreach() {
        thl[] = thvc*H[] + T_ref*(1 - H[]);
        qt[]  = qtc*H[]  + qt0*(1 - H[]);
        u.y[] = wc*H[];
        u.x[] = 0.001*noise();
      }
      TOLERANCE = 1e-4;
    }

    The temperature and humidity are subject to diffusion and a in-cloud forcing.

    event tracer_diffusion(i++) {
      diffusion (thl, dt, mu);
      diffusion (qt, dt, mu);
      foreach() {
        thl[] += H[]*dt * (thvc - thl[])/tau;
        qt[]  += H[]*dt * (qtc  - qt[] )/tau;
      }
    }

    The effects of gravity are including via the prognostic buoyancy field. Further, the velocity in the cloud is subject to a forcing.

    event acceleration (i++) {
      scalar b[];
      foreach()
        b[] = B_CENTERED;
      boundary ({b, H});
      foreach_face(y)
        av.y[] = ((b[] + b[0,-1])/2. +
    	      (wc - (u.y[] + u.y[0,-1])/2.)/tau * (H[] + H[0,-1])/2.);
    }

    Grid Adaptation

    In order to focus the computational resources, we adaptively refine and coarsen the grid based on the discrete representation of the temperature, humidity and velocity-component fields.

    event adapt (i++) {
      adapt_wavelet ({thl, qt, u}, (double[]){6./10., 0.005/10., wc/10., wc/10., wc/10.}, 9);
    }

    Results

    We output a movie:

    Looks OK

    movie generation code:

    #include "view.h"
    
    void cloud (double cmap[NCMAP][3]) {
      for (int i = 0; i < NCMAP; i++) {
        cmap[i][0]= 0.2 + 0.8 * (double)i/(double)NCMAP;
        cmap[i][1]= 0.2 + 0.8 * (double)i/(double)NCMAP;
        cmap[i][2]= 0.9 + 0.1 * (double)i/(double)NCMAP; 
      }
    }
    
    event movie (t += 0.5) {
      
      scalar b[], ql[];
      foreach() {
        b[] = B_CENTERED;
        ql[] = QC;
      }
      boundary ({ql, b});
      
      view (fov = 40, width = 900, height = 900);
    
      squares ("b", min = -0.05, max = 0.05, linear = true);  
      mirror ({-1,0}) {
        cells();
      }
      translate (y = -L0) {
        squares ("u.y", min = -wc, max = wc, linear = true);  
        mirror ({-1,0}) 
          squares ("ql", map = cloud, min = 0, max = 0.002, linear = true);
      }
      draw_string ("Buoyancy", 2);
      draw_string ("Grid Cells", 1);
      draw_string ("Vertical velocity", 3);
      draw_string ("Liquid water", 0);
      save ("sss.mp4");
    }
    
    event stop (t = 120);

    Reference:

    Nair, V., T. Heus, and M. van Reeuwijk, Dynamics of Subsiding Shells In Actively Growing Clouds With Vertical Updrafts. J. Atmos. Sci., https://doi.org/10.1175/JAS-D-19-0018.1