Embedded Files
OVERVIEW
OVERVIEW
In this chapter we'll take a look at a few key processes in warm-cloud microphysics: (i) aerosol activation-condensation-evaporation, (ii) collision-coalescence. The activation-condensation-evaporation solves a number of coupled ODEs. At this point we're familiar with ODE solvers so we'll just use scipy's solve_ivp. Collision-coalescence looks very similar to "chemical reactions" but the calculation is performed in a nested loop over all size bins.
Here we use a "sectional approach", i.e. we track aerosols and droplets in a number of bins. The math is simpler because it's based on first principles, but is computationally expensive. As a result, most cloud microphysics schemes used in 3-D models don't use sectional approach. Instead, "bulk" schemes are commonly used in 3-D models because of the computational efficiency but bulk schemes have built-in assumptions. For instance, 2-moment bulk schemes prognose mass and number of hydrometeors but not size (at least not explicitly).
The full Jupyter Notebook is hosted on Github: YAMCHA_tutorial_ch12.
1. Droplet Growth Rate
1. Droplet Growth Rate
Firstly let's calculate the ambient supersaturation (S), which is the ratio between specific humidity (qv) and saturation water vapor concentration (qsat) minus 1.
S = qv / qsat - 1
Now let's look at the first differential equation describes the growth rate of the droplet diameter. D_wet(k) is the diameter of one droplet in the k-th bin, S is the ambient supersaturation, S_eq(k) is the equilibrium supersaturation of the k-th bin.
dD_wet(k)/dt = 4G(k)*(S - S_eq(k))/D_wet(k)
G(k) is a condensation coneffidient of the k-th bin, which consists of two terms: vapor diffusion and heat transfer. rho_w is water density, Dv_eff is the effective water vapor diffusivity in air corrected for non-continuum effect, es is the saturation vapor pressure of water, Rv is ideal gas constant divided by molecular weight (of water), Lv is the latent heat for water condensation, Ka_eff is the effective thermal conductivity of air corrected for non-continuum effect.
G(k) = (vapor_diffusion + heat_transfer)^-1
vapor_diffusion = rho_w*Rv*T/Dv_eff/es
heat_transfer = rho_w*Lv*(Lv/Rv/T-1)/Ka_eff/T)
Let's break down. S - S_eq(k) describes the thermodynamic driving force. S > S_eq(k) means the surrounding environment containts more water vapor than the droplet can support at equilibrium, i.e. condensation will occur and the droplet will grow. S < S_eq(k) means the surrounding environment is "too dry" relative to the droplet, i.e. evaporation will occur and the droplet will shrink.
G(k) describes the kinetics, i.e. how fast the droplet responds to the thermodynamics. G(k) has two terms describing the vapor diffusion rate and heat transfer rate. When a droplet is growing, water vapor must diffuse from the surrounding environment to the droplet surface, while the latent heat released from condensation must be conducted away. Both can affect how fst droplet can grow.
2. κ-Köhler Theory
2. κ-Köhler Theory
qsat is the saturation vapor concentration over flat surface, which can be approximated by a nice little function of temperature. Saturation vapor pressure over curved a surface such as a droplet can be very different. κ-Köhler Theory describes the supersaturation over droplets at equilibrium S_eq(k), i.e. related to the saturation vapor pressure (concentration) over droplets. S_eq(k) is affected by two factors:
S_eq(k) = Kelvin * Solute - 1
Kelvin effect: vapor pressure increase over curved surfaces, which depends on surface tention.
Kelvin = exp(4*SurfaceTension*Mw/R/T/rho_w/D_wet(k))
Solute effect (essentially Raoult's law): dissolved solutes lowers the vapor pressure, making condensation easier. κ is the hygroscopicity of solute.
Solute = (D_wet(k)^3 - D_dry(k)^3)/(D_wet(k)^3 - (1-κ)*D_dry(k)^3)
Hygroscopicity (κ) depends on the composition of the aerosol. Higher κ means stronger tendency to take up water (e.g. sea salt). Figure on the right demonstrates these effects. Kelvin effect greatly limits the growth of smaller droplets (need much higher super saturation). Solute effect can promotes growth by reducing supersaturation.
Now let's code it up! Fairly straightfoward. You can use it to make that lil plot showing how dry diameter and hygroscopicity affect supersaturation.
3. Vapor-Liquid Conversion: Mass Conservation
3. Vapor-Liquid Conversion: Mass Conservation
Say, the mass of a droplet in the k-th bin is m(k). The surface area of the droplet is pi*D_wet(k)^2. With a droplet diameter growth rate of dD_wet(k)/dt, the radius change is just half of that. Consider a very small increase in radius (delta_r), then the change of volume is just surface_area * delta_r. Therefore, the volume change rate of this droplet is pi*D_wet(k)^2 * 0.5* dD_wet(k)/dt. With density, the change of mass m(k) is given by:
dm(k)/dt = rho_w*pi*D_wet(k)^2 * 0.5* dD_wet(k)/dt
This is the change of liquid water mass of one single droplet in k-th bin. The droplet number concentration is N(k). Sum up the liquid water mass change in all bins, that's how much liquid water is condensed, which is also how much water is removed from the vapor-phase:
dqv/dt = -sum(N(k)*dm(k)/dt)/rho_air
These 2 differential equations ensure mass conservation for total water (vapor + liquid).
4. Temperature
4. Temperature
Last but not least, the temperature tendency is driven by both adiabatic cooling (1st term on the RHS) and latent heat from condensation (2nd term on the RHS, depends on dqv/dt). w is the updraft speed, Cp is the specific heat capacity of air.
dT/dt = -g/Cp*w - Lv/Cp*dqv/dt
g/Cp is just the dry adiabatic lapse rate.
Note that entrainment is ignored for now. In reality, entrainment drives the mix between the air parcel and the environment.
5. Pressure and Altitude
5. Pressure and Altitude
z is the altitude of the parcel. With an updraft speed w, the change of height dz/dt is just:
dz/dt = w
Assume hydrostatic balance for the air parcel, i.e. pressure gradient (dp/dz) is balanced by gravity (rho*g):
dp/dz = -rho_air*g
dp/dt = -rho_air*g*dz/dt = -rho_air*g*w
These two additional differential equations complete the simple parcel model! Again we don't have entrainment at this point.
6. Derivatives for the Solver
6. Derivatives for the Solver
Now the fun part! We have 6 coupled ODEs. From previous chapters we have learned how to use scipy's ODE solver. We can absolute use that to sove these ODEs.
dD_wet(k)/dt, dm(k)/dt, dqv/dt, dT/dt, dp/dt, dz/dt
Before going too far, let's think about this. Dwet(k) and m(k) are arrays, i.e. for each bin k, we'll have a wet diameter D_wet and a mass m. If you have 10 bins (often not enough), then you'll have 10+10+4 coupled ODEs. Good thing is, we don't necessarily need to prognose per-droplet mass m(k), because we can diagnose that as long as we have wet diameter (assume droplets are perfect spheres). So the problem is simplified.
dD_wet(k)/dt, dqv/dt, dT/dt, dp/dt, dz/dt
Remember if we want to solve the ODE, we need to assemble the RHS derivative function. For our chemistry problems, we use a parser to generate that automatically. But for the parcel model, it's simple enough, so we'll just manually code that:
Here all state variables (e.g. D, qv, T, p, z) are packed into one common array y. This is easier for the solver to handle but not so much for us humans. Therefore, at the very beginning of this function, we unpack y into local variables Dwet_um, qv, T, p, and z. The rest of this function should be straightforward.
A few important notes:
This solves aerosol activation, condensation, and evaporation too. None of these affect the aerosol number in each bin N(k), i.e. aerosol/droplet number is conserved. This is not the case for other processes, e.g. collision-coalescence (will cover this later).
We track the growth/shrink of each individual aerosol/droplet bin. This is one common approach, often called "sectional" or "bin-based" approach, as appose to "bulk" or "modal" approaches. Generally speaking, "sectional" approaches are more accurate but are computationally more expensive.
7. Solve the ODEs!
7. Solve the ODEs!
Before calling the solver we need to specify a few more things, including the initial size distribution. We'll use a function that we defined in Chapter 8.
# --- integrate dN/dDp in certain rangesdef integ_dN_dDp_range(in_dN_dDp, in_Dp, in_Dp_lower, in_Dp_upper): Dp_lower_ind = abs(Dp_um-in_Dp_lower).argmin() Dp_upper_ind = abs(Dp_um-in_Dp_upper).argmin() return cumulative_trapezoid(in_dN_dDp[Dp_lower_ind:Dp_upper_ind], in_Dp[Dp_lower_ind:Dp_upper_ind])[-1]
# --- calcualte bin bounds from bin centersdef get_bin_bounds(in_Dp_um): out_Dp_bounds_um = np.zeros(1+len(in_Dp_um)) for ibin in range(len(in_Dp_um)): if ibin==0: out_Dp_bounds_um[1] = 0.5*(in_Dp_um[ibin]+in_Dp_um[ibin+1]) out_Dp_bounds_um[0] = in_Dp_um[ibin]-(out_Dp_bounds_um[1]-in_Dp_um[ibin]) elif ibin==-1+len(in_Dp_um): out_Dp_bounds_um[ibin] = 0.5*(in_Dp_um[ibin]+in_Dp_um[ibin-1]) out_Dp_bounds_um[ibin+1] = in_Dp_um[ibin]+(in_Dp_um[ibin]-out_Dp_bounds_um[ibin]) else: out_Dp_bounds_um[ibin] = 0.5*(in_Dp_um[ibin]+in_Dp_um[ibin-1]) return out_Dp_bounds_um
Let's define the initial size distribution as a lognormally distributed function, 51 bins. Median diameter is 1.5 micron and sigma is 2.
rho_aerosol_kg_m3 = 1800.init_N_cm3 = lognormal_dist_dN_dDp(D_um, in_Dp_mean=1.5, in_sigma=2., in_N0=0.5)init_dN_dlogDp_cm3 = init_N_cm3/np.diff(np.log10(D_bounds_um))init_dV_dlogDp_um3_cm3 = init_N_cm3*(np.pi/6.)*(D_um**3.)/np.diff(np.log10(D_bounds_um))
With the RHS derivative function, we can now call the solver!
# --- set timetime_s = np.arange(0,300,dt)dt = 10
# --- call the solversol = solve_ivp(parcel_rhs, t_span=(0, time_s[-1]), y0=np.concatenate((D_um,[qv0,T0,p0,0])), t_eval=time_s, args=(D_um, init_N_cm3, 0, 1.1), method="BDF",rtol=1e-6,atol=1e-10)
# --- get the new sizes!D_wet_um = sol.y[:-4,:]
# --- liquid water contentLWC_g_m3 = [np.sum(init_N_cm3*(np.pi/6.)*(D_wet_um[:,i]**3.)) * 1e-6 for i in range(len(time_s))]
Here y0 is the initial condition of the state variable y, as you can see we pack D_um, qv0, T0, p0, and z0 (set to zero here) into one array. 4 args are passed into the derivative function: dry diameter (D_um), number concentration (init_N_cm3), updraft speed (set to zero here), as well as kappa (set to 1.1 here). At the end of the simulation, we get the the new diameters, D_wet_um. Do you know what the last 4 columns in sol.y are?
In this test we run the model for 300 seconds with time step being 10 seconds. The plot on the right shows the simulation results at 95% and 100.1% relative humidity.
8. Collision-Coalescence
8. Collision-Coalescence
Now we'll switch gear and look at another very important process. From the previous equations you have probably realized that the condensational growth is less efficient for large droplets. Indeed, with only condensation, droplet growth beyond 10-20 μm is inefficient.
Above ~20 μm, collision-coalescence becomes increasingly important. Collision-coalescence describes smaller cloud droplets bump into each other and merge to grow large droplets (similar to aerosol coagulation). A main driver is gravitational coalescence, i.e. larger droplets fall faster than smaller droplets, so smaller droplets are collected by large droplets.
The prognostic equation for collision-coalescence is given in the figure. dN(i)/dt has two terms: the gain of droplets in bin i due to collisions between bins j and k (source), also the loss of droplets in bin i due to collision with bin j (sink). Don't these look at chemical reaction rates? Thats exactly it! For instance, the collision rate between bins j and k is:
K(j,k)*N(j)*N(k)
Note when implementing the "source" term, there is a factor of 1/2. This is because of the nested loops, i.e. the collision between bins j and k will be counted twice. This factor of 1/2 is to avoid double counting.
Here the "rate coefficient" K is called coalescence kernel, which depends on collection efficiency E, cross section, and relative terminal velocities.
Collection efficiency: depending a lot on relative sizes. Here we use the look-up table in Hall (1980), in which collection efficiency is given as a function of the ratio between big and small droplets.
Terminal velocity: the highest speed when an object falls! Also depends strongly on size. Note that if you use the gravitational setting we discussed before, it'll drastically overestimate for droplets larger than ~100 microns. We use the parameterization in Beard (1977).
Here's my collision-coalescence implementation (full version in the jupyter notebook):
IMPORTANT NOTES:
This uses a simple forward Eulerian solver so the time step shouldn't be too big.
Collision-coalescence does alter number concentration! This is different from condensation/evaporation.
This approach does not alter size distribution, i.e. new bin produced from collision is mapped into the nearest 2 bins in the original size bins (see the codes highlighted in red).
In the Jupyter notebook, we have a test case which is configured based on Bott (1998), and the result is given in the figure below. Note Bott (1988) used a different parameterization for size distribution. Also here we use diameter instead of radius. But the results look pretty similar to Fig 3(b) in Bott (1998), i.e. it successfully produced a drizzle mode that extends to >1000 micron.
When large droplets are falling, they may break up and produce smaller droplets. This is increasingly important for droplets >1000 microns (>1 mm).
9. Final Notes
9. Final Notes
In this chapter, we demonstrate activation-condensation-evaporation and collision-coalescence separately. Activation-condensation-evaporation is solved on moving-bins (also called Lagrangian-bin), while collision-coalescence is solved on fixed bins (mapped to the original bins). These are perhaps the simplest way to implement these processes. Moving-bin approach is pretty accurate. Fixed bins, however, may broaden/widen modes which is a numerical artifact.
In the activation-condensation-evaporation scheme, we prognose size but not mass or number, because number is not affected and liquid water content can ge diagnosed with wet size. Water condensation/evaporation does not alter aerosol mass either, so we don't need to prognose aerosol masses here either.
Collision-coalescence is a lot more complicated, because it alters number, mass of liquid water, also mass of aerosols. In this simple demo, we prognose only numbers but not mass. If you take another look at the codes, we're pretty close to prognosing liquid water mass already. You can also modify the codes to carry aerosol mass in each bin. Both scale with number. Needless to say, mass of water and aerosols must be conserved (important sanity checks). For simplicity we won't dive into this.
These are probably the most important processes in warm-cloud microphysics! A few processes not covered here are:
Gravitational setting: easy to implement. You know how to calculate gravitational setting velocity/terminal velocities
Large droplets breaking up while falling: starts to become important for droplets >~1mm
Entrainment: easy to implement as long as you have the entrainment rate, but getting the entrainment rate can be tricky
Aqueous-phase chemistry: we discussed explicit phase-transfer and aqueous-phase chemistry in Chapter 7 already
All processes discussed here are driven by "mean atmospheric motions". In reality, turbulence can drive cloud processes as well
These are just physical processes. To use these in an application, you need a host model framework. There are several ways to do this. The simplest is to include these in a 0-D parcel model, you can then force the parcel to ascent/descent or force a change in supersaturation. You can also implement these in a 1-D framework (e.g. Chapter 11), obviously this will be more complicated, especially if you want to implement a vertical velocity. The 1-D framework demonstrated in this tutorial is "old fashioned": strictly speaking, 1-D models don't deal with vertical velocity which is a result of horizontal wind convergence/divergence which are beyond the scope of 1-D models. But you can write an advection solver (w*dC/dz) and add it on top of the diffusion solver (d/dt(K*dC/dz)), i.e. turn it into an advection-diffusion solver (similar to the Burger's equation). You can then prescribe a vertical velocity profile (w).
10. Aerosol Mixing State?
10. Aerosol Mixing State?
Up to this point, all aerosol related modules in this tutorial assumes aerosols are externally mixed, i.e. each individual aerosol particle consists of a single, pure chemical species or type. In fact, most (maybe all) 3-D models assume aerosols are externally mixed. This is often a gross oversimplification. If you look at actual images of aerosol particles, you'll often see organics and sulfate coated on black carbon "cores", or dust co-existing with sea spray. But when these black carbon or dust or sea spray aerosols are emitted, they are "pure". How do they become internally mixed with other stuff?
Several ways can produce internally mixed aerosols, such as condensation of organic vapor or sulfuric acid onto black carbon "cores" (chemistry), or coagulation of aerosols. Collision-coalescence may also produce internally mixed aerosols! For instance, when sea spray aerosols are activated into droplets, they can collect with dust or tiny droplets produced from dust (dust are also hygroscopic, just less so compared to other types). Collision-coalescence involving externally mixed sea spray and dust can produce droplets containing both sea spray and dust; after water is evaporated, aerosols containing sea spray and dust would stay in the particle-phase, producing internally mixed sea spray-dust aerosols!
A modified version of the collision-coalescence scheme is used to calculate the production of internally mixed sea spray-dust aerosols and the preliminary results are shown on the right. This scheme is a lot more complicated than the demo shown here, because it has to track not only number, but also masses of individual particles before collision. But this demonstrates how cloud processing can affect the composition of aerosols!
Google Sites
Report abuse