Siyuan Wang | Research - Ch 10. Vertical Diffusion

Ch 10. Vertical Diffusion

OVERVIEW

- In all previous chapters the model is zero-dimensional, a.k.a. "box model". Everything is homogeneously mixed immediately inside the box. You may add boundary condition but it's static. In this chapter we'll have more than one box!
- We'll have a number of boxes lining up in one direction, and we'll model the exchange between these boxes. This modeling concept, one-dimensional (1-D) modeling, is quite useful, for instance, you can build a 1-D model to study the development of atmospheric boundary layer and the impact on air pollution, you can also build a 1-D model to mimic the multi-layer structure of aerosols and its impact on aerosol chemistry. In fact, many three-dimentional (3-D) chemistry transport models are "column-based", i.e. physics are solved in a number of columns, and the exchange between those columns are taken care of by horizontal advection.
- One process that affects the exchange between the boxes, or "layers", in a 1-D model, is diffusion. The diffusion equation involves 2nd order derivatives, therefore cannot be directly solved by ODE solvers as chemistry. We will learn two approaches to solve the diffusion equation. Note that this chapter does not cover advection.
- We'll also go over how to use SciPy's built-in ODE solver, which has been used extensively in the chemistry solver and in prevoius Chapters.

1. The Diffusion Equation: Fundamentals

The 1-D model we will be building has a number of layers stacking vertically on top of each other. The diagram on the right shows what these layers may look like. One important thing is that, in 1-D modeling, some variables are represented at the center of the layers, and some are at the interfaces between the layers, or at the top/bottom of the model. Obviously the number of interfaces, including the top and bottom interfaces, is the number of layers plus one.
The prognostic equation of 1-D diffusion is:

dc/dt = K*d2c/dz2

where c is the tracer concentration which is a function of time (t) and height (z). K is the diffusivity, and its SI unit is m2/s. For now let's assume K is constant.

We'll solve this using finite difference method, in a uniform z grid with spacing dz (layer thickness). The continuous function is represented by a discretized function:

f[zi] = c[i]

where zi is the height of i-th layer, and c[i] is just the tracer concentration at this i-th layer. It's Taylor serios at zi is then:

f[zi+dz] = f[zi] + sum( dz^j/j! * f^(j) )

Please excuse my notation but here f^(j) is the j-th order derivate of f, so f^(1) is just f', f^(2) is f'', etc.

To get the 2nd order derivative, run the first few terms of the Taylor expansion forward and backward around zi:

f[zi+dz] = f[zi] + dz*f'(zi) + (dz^2/2!)*f''[zi] + (dz^3/3!)*f'''[zi] + ...

f[zi-dz] = f[zi] - dz*f'(zi) + (dz^2/2!)*f''[zi] - (dz^3/3!)*f'''[zi] + ...

Add these together, we get:

f[zi+dz] + f[zi-dz] = 2*f[zi] + (dz^2)*f''(zi) + ...

f''(zi) = (f[zi+dz] - 2*f[zi] + f[zi-dz]) / (dz^2)

d2c/dz2 = (c[i+1] - 2*c[i] + c[i-1]) / (dz^2)

High order term dropped. This is the approximation of the 2nd order derivative that we can use in our numerical schemes!

dc/dt = K * (c[i+1] - 2*c[i] + c[i-1]) / (dz^2)

Previously we assumed K is a constant, which is fine if, say, you want to build a multi-layer aerosol model. But if you want to simulate the vertical diffusion of compounds in the atmosphere, K certainly varies with height. Let's call this Kz, i.e. height-dependent diffusivity.
Rearrange the prognostic equation of 1-D diffusion we get:

dc/dt = (d/dz)*(Kz*dc/dz)

This term, Kz*dc/dz is actually teh flux of the tracer. Let's call it F.

Flux is normally solved at the layer interface rather than the layer center, so it's easier to make sure the mass is conserved. Therefore:

dc[i]/dt = dF/dz = (F[i] - F[i+1])/dz

dc[i]/dt = (Kz[i]*(c[i-1]-c[i])/dz - Kz[i+1]*(c[i]-c[i+1])/dz)/dz

where F[i] is the flux across the lower interface of i-th layer, which is the same flux acrossing the upper interface of i-1 th layer. Similarly, F[i+1] is the flux across the lower interface of i+1 th layer, which is the same flux across the upper interface of i th layer. See the plot on the right!

Note that we also need to take care of the boundary condition, i.e. what we'll do with the upper most and lower most interfaces/layers. There are several different types of boundary conditions, e.g. you can have the so-called "periodic boundary" in which case whatever leaves one side of the model column comes back in from the opposite side, you can prescribe certain concentrations (fixed or varying) in the top/bottom most layers, or prescribe certain fluxes across the top/bottom most interfaces. The latter has a special case, with zero fluxes across the top/bottom most interfaces. Obviously if we have zero fluxes across the top/bottom most interfaces, all species are "trapped" inside the column. With no chemical reactions, the total mass will be conserved. This is what we'll do today.

2. Solve as Ordinary Differential Equations (ODEs)

Grab the all-in-one Jupyter Notebook (YAMCHA_tutorial_ch10) from GitHub.
Run the first Cell and load all the libraries. Navigate to the 2nd cell. At this point you're familiar with SciPy's solve_ivp function already, which we have been using extensively to solve chemistry. The diffusion equation cannot be directly solved using ODE solvers, but we have approximated the 2nd order derivative term so it can be solved as ODE! We just need to write a function that calculates the derivatives:

def dcodt_vertical_diffusion(t,conc):

dcodt_out = np.zeros(len(z_m))

for i in range(0,nz):

if i==0: # bottom of the model

dcodt_out[i] = (1./dz_m/dz_m)*(-kz_m2_s[i+1]*(conc[i]-conc[i+1]))

elif i==-1+nz: # top of the model

dcodt_out[i] = (1./dz_m/dz_m)*(kz_m2_s[i]*(conc[i-1]-conc[i]))

else: # all other layers

dcodt_out[i] = (1./dz_m/dz_m)*(kz_m2_s[i]*(conc[i-1]-conc[i])-kz_m2_s[i+1]*(conc[i]-conc[i+1]))

return dcodt_out

The line highlighted in red is the generic implementation of the discretized form of the prognostic equation. The other two if/elif statements take care of the boundary condition. In this case, in the bottom layer, we only have the flux going up to the next layer, and no flux coming in from the bottom interface.

3. Set up the Height Grid

Now let's setup the height grid, and create the layer center height (z_m) and layer interface height (z_interface_m). The layer thickness, dz_m, is set to 100 meters here, and the model top goes up to 3000 meters. This configuration has 31 vertical layers.

# --- height grid: layer center

dz_m = 100

z_m = np.arange(dz_m/2, 3000.+dz_m, dz_m)

nz = len(z_m)

# --- height grid: layer interface

z_interface_m = np.zeros(len(z_m)+1)

for i in range(len(z_interface_m)):

if i==0: z_interface_m[i] = 0

else: z_interface_m[i] = 2.*(z_m[i-1]-z_interface_m[i-1])+z_interface_m[i-1]

4. Set up the Height-Dependent Diffusivity

This is a fairly complicated concept but is covered in every meteorology text book. I'm not going to cover too much details here, we'll just generate a reasonable vertical profile of diffusivity, Kz, and use that to demonstrate the 1-D diffusion.
The Kz profile I got is shown on the right. As you can see it starts at zero at the surface and peaks at ~120 m2/s at ~500 m, then drops with increasing altitude to about 1500 m. With this profile, we'd expect considerable mixing in the lowest ~1500 m (we'll call it mixing layer) and significantly reduced mixing above ~1500 m.

5. Run the Model!

...which are fairly straightforward. The only thing I should mention is the initial condition:

# --- initial condition

conc[0,20] =1e+7

conc[0,5] =1e+7

Basically, the initial profile of the tracer has two peaks, one in the 5th layer (within the mixing layer) and the other in the 20th layer (above the mixing layer), both at 1e+7 molec/cm3.

We're running the model for 0.3 days and saving outputs every 10 seconds. With only one tracer and only 31 layers, it took me about 10 seconds to complete the simulation. Not bad!
The following plots show the results. The plot on the left shows a few vertical profiles: the blue line is the initial profile with two sharp peaks, and all other colorful lines are just profiles at different times. We can see that the peak within the mixing layer (<1500 m) undergoes major diffusion, eventually well mixed throughout the mixing layer. The peak emitted above the mixing layer (~2000 m), on the other hand, undergoes very limited vertical diffusion. The curtain plot on the right shows the same thing!

Mass conservation! Very important. Can you spot where in the Jupyter Notebook do I perform this important sanity check?

6. Solving Diffusion using ODE Solver: ODE Fundamentals

Here we solve the diffusion equations by discretization and use an ODE solver. Might as well take the chance to say a few more words about these ODE solvers! The SciPy documentation is pretty straightforward. The most important thing you need to provide is the function that contains the derivatives. Then you need to take care a few more things in the solver interface. For instance:

conc[i,:] = solve_ivp(dcodt_vertical_diffusion,

tout_s[0:2],

conc[i-1,:],

t_eval=[output_frequency_seconds],

method='BDF',

rtol=1e-5,atol=1e-5).y[:,0]

This is the line where the the solver is called. A few important things here:

Derivative function: it's just this guy dcodt_vertical_diffusion that we defined previously.

Integration time range: here given as tout_s[0:2], but since we're grabing the first two time steps in the list only, it's equivalent to [0, output_frequency_seconds. You can provide the entire time range here (tout_s) but we don't do it because the condition/kinetics may well change between timesteps.

Initial condition: it's just the results from the previous step (conc[i-1,:]). As a result, at the first timestamp (t = 0), the solver is not called.

Method: very important. Here because I copy-pasted this from our previous chemistry models, in which the chemical systems are usually extremely stiff. For stiff problems, 'BDF' or 'Radau' is better, 'LSODA' as a general one-fit-all solver also works but it's a bit slow. For non-stiff problems, start with 'RK23' or 'RK45'.

Tolerance: these numerical approaches are not designed to get the "right" answers, but rather, they look for the answer that is "just right". You specify some tolerance level, and as long as the results are within the tolerance level, the solver calls it good and moves on. Here tolerance is controled by both relative and absolute (rtol and atol), here are both set to 1e-5. If you run into numerical issues or the solver is too slow, you can always relax the tolerance.

The nice thing about using the ODE solver for diffusion is that you can couple it with chemistry! Coupling chemistry directly with the diffusion equations is straightforward also accurate but may not be the best idea, especially if you have multiple species (hundreds or thousands) or many vertical layers. The computational cost ramp up quickly.
Let's do a quick test: by default, we'll stick to the 'BDF' method, and let's increase the vertical resolution of the model and change dz_m from 100 m to 30 m. Now the number of vertical layers is increased by 3x and we got 101 vertical layers! Run the model again. The time needed to complete the simulation increased dramatically: from ~10 seconds to 54 seconds! If switch to 'RK23' method, with 100 m resolution it takes 2.6 seconds to complete! Even with 30 m resolution it only takes 8 seconds! Much faster than the 'BDF' method! This is because the diffusion problem, with the 2nd-order derivatives discretized, isn't particularly stiff, therefore an explicit method like the Runge-Kutta (RK) series is usually faster.
Next, we'll explore other approaches to solve the diffusion equations by using the numerical scheme we derived in the beginning!

7. Solving Diffusion using Explicit Euler Scheme

Previously we got the the discretized prognostic equation.

dc[i]/dt = (Kz[i]*(c[i-1]-c[i]) - Kz[i+1]*(c[i]-c[i+1]))/dz

and if we take a small time step, the left hand side is just:

dc[i]/dt = (c[n,i] - c[n-1,i]) / (t[n]-t[n-1])

where c[n,i] is the concentration at n-th time step at the i-th layer...

This is very easy to implement in the model. Navigate to Cell #3, and the following snippet does just that:

for n in range(1,n_tstep):

for i in range(0,nz):

if i==0: # bottom of the model

conc[n,i] = conc[n-1,i] + dt/dz_m/dz_m*(-kz_m2_s[i+1]*(conc[n-1,i]-conc[n-1,i+1]))

# add emissions at the bottom!

# conc[n,i] = conc[n,i] + emis_flux_molec_cm2_s[i]/(dz_m*100.)*dt

elif i==-1+nz: # top of the model

conc[n,i] = conc[n-1,i] + dt/dz_m/dz_m*(kz_m2_s[i]*(conc[n-1,i-1]-conc[n-1,i]))

else: # all other layers

conc[n,i] = conc[n-1,i] + dt/dz_m/dz_m*(kz_m2_s[i]*(conc[n-1,i-1]-conc[n-1,i])-kz_m2_s[i+1]*(conc[n-1,i]-conc[n-1,i+1]))

See the line highlighted in red. The other two if/elif statements take care of the boundary condition.

Run the model! With dz = 100 m, it takes 0.17 seconds to complete, and the results are nearly identical to that using the ODE solver! Even with dz = 30 m, it takes 0.6 seconds to complete, much faster than the ODE-based scheme even with 'RK23'.

8. Explicit Scheme: Stability Concerns

However, this explicit scheme has stability concerns! For instance, what happens if you further decrease dz to 20 m, and keep everything unchanged (especially dt = 10 seconds)? It took <1 second to finish but produced some totally bananas results. See the figure on the right.
In this case, dz = 10 s and dz = 20 m, and the model becomes instable and leads to instable and oscillating results. This is an important consideration of explicit schemes like this. In layman's terms, the diffusion "speed" resulted from the flux is too fast for the grid spacing, dz!
An quick and easy fix is to decrease the time step, dz. For instance, previously the time dz is 10 s; if you just decrease dz to 5 s or 1 s, the model still runs very fast but will be able to produce stable results.
If you run sophisticated 3-D models (e.g. WRF), you may have run into issues with error messages telling you that there are CFL violations! Essentially it's the same idea. If you run the model at higher spatial resolution, you need to decrease your time step accordingly. This is why high resolution models are exponentially more expensive than low resolution models.

9. Add Emissions as Lower Boundary Condition

We mentioned before that one way to provide boundary condition for the 1-D model is to specify the fluxes coming in/out of the model at the top/bottom interfaces. Here I'm going to demonstrate how this can be done.
You may notice a line where I define a 1-D array called emis_flux_molec_cm2_s, which has a constant value of 1e+7, with the unit of molec/cm2/s. Then in the big loop, a line was commented out where the lower boundary condition is defined. Decomment that line:

conc[n,i] = conc[n,i] + emis_flux_molec_cm2_s[i]/(dz_m*100.)*dt

As you can see, essentially we convert the emission flux (molec/cm2/s) to emission rate (molec/cm3/s) by dividing the layer thickness (dz_m, need to convert the unit from m to cm), which is then multiplied by the time step (dt) and added to the concentration after multiplying. Fairly straightforward!
Before running the model, let's revert the layer thickness dz_m back to 100 m. Then run the model and the curtain plot with emissions is shown here! You can see that at the bottom of the model the tracer concentration starts to pick up and accumulate inside the mixing layer. If you run some mass conservation test you can retro calculate how much emission you have injected!

10. Summary

In this chapter, we learned how to use SciPy's built-in ODE solver, atlhough you have been using it in previous chapters a lot already!
Then we learned some fundamentals of 1-D diffusion, and used two different approaches to solve the diffusion equations: an ODE solver, and an explicit Euler scheme. The ODE solver is straightforward but quite slow. The explicit scheme is a lot faster but you need to worry about stability.
1-D models are not only a powerful tool in atmospheric science, but also a cornerstone in 3-D models. With the same diffusion solver, you can model other processes as well! For instance, the prognostic equations for heat and moisture are:

dT/dt = K*d2T/dz2 + S_heat

dq/dt = K*d2q/dz2 - Cd

dqc/dt = K*d2qc/dz2 + Cd

where T is temperature (actually potential temperature is often used), K is the turbulent diffusivity for heat (also used for moisture and other tracers), S_heat is additional heat source (e.g. sensible heat, or latent heat released from condensaton), q is the specific humidity, and qc is the condensed-phase water concentrations (e.g. cloud), and Cd is the condensation rate. The moisture processes (e.g. how water vapor and condensed-phase water evolve) are usually solved in sophisticated microphysics models. But you can parameterize simple condensation rate and code up your own simple microphysics model! Couple that with the heat equation and, bang, you got yourself a 1-D model that is comprehensive enough so it doesn't require temperature and specific humidity input as boundary conditions!

Google Sites

Report abuse