Lab 5-3: Winds and Orographic Precipitation#
Our goal here is to combine what we learned in Lab 2 regarding atmospheric thermodyanamics to better understand rain rates along a mountain. This lab draws on the Precipitation lecture notes, as well as on a paper written by Gerard Roe in 2005.
import numpy as np
import pandas as pd
import scipy.stats as stats
import matplotlib.pyplot as plt
%matplotlib inline
Note that wind calculations will be added in a future iteration, for now wind information will be given.
Saturated vapor pressure is a function of air temperature. In lab 2 (and in the Shuttleworth book), we used the formula, as shown here, to calculate saturated vapor pressure in units of kPa with an input of temperature in degrees C. $\( e_{sat} = 0.6108*exp(\frac{17.27T}{237.3+T}) \)$
In Gerard Roe’s paper, he describes this equation with the Clausius-Clapyron relationship, which he writes as $\( e_{sat} = 6.112*exp(\frac{aT}{b+T}) \)\( where here \)e_{sat}\( is in millibars or hPa, a = 17.67, and b=243.5\)^\circ$C.
Now, it’s always somewhat frustrating when multiple reputable sources have different equations for the same thing. It’s a good idea to pause and check to see if they’re really quite close to the same thing. One way to do this is a quick plot.
T_curve=list(range(-40,40))
T_curve=np.asarray(T_curve)
# We turn it into a numpy array so that we can use numpy to calculate esat for the entire array
esat1 = 0.6108*np.exp((17.27*T_curve)/(237.3+T_curve))
# this formula outputs esat in kPa, if we multiply by 10, we get hPa, easier for comparison
esat1 = esat1*10
a=17.67
b=243.5
esat2 = 6.112*np.exp((a*T_curve)/(b+T_curve))
# Create a new figure.
plt.figure()
# Use the plot() function to plot the year on the x-axis, peak flow values on
# the y-axis with an open circle representing each peak flow value.
plt.plot(T_curve, # our x value
esat1 # our y value
)
plt.plot(T_curve, # our x value
esat2 # our y value
)
# Label the axes and title.
plt.xlabel('Temperature (deg C)')
plt.ylabel('Saturated Vapor Pressure (kPa)')
plt.title('(Fig. 1) Saturated Vapor Pressure Curve');
Try commenting out one line in the above code so that only one plots. These are virtually indistinguishable. They are different numerical approximations of the same thing. Wikipedia has an interesting review of a number of other formulas and their accuracy at different air temperatures, if you are interested. In many cases, errors arise at high temperatures (e.g., for applications in pressure cooking or boilers). For the purposes of our class, either of the formulas here is fine.
Pause and think#
The plot above is very important to our considerations of precipitation, storms, floods, and flood probabilities. The first step for precipitation to form is for the atmosphere to carry enough moisture (water vapor) to reach the saturated vapor pressure long enough for condensation to form a lot of droplets, which fall as precipitation. As we discussed in class, the heaviest rain events in the Pacific Northwest are from atmospheric rivers, which are warmer than average storms and carry a lot of moisture from the tropics in a narrow band.
The plot above should also give you pause as you think about warming temperatures. Within typical atmospheric temperatures, the saturated vapor pressure increases rapidly with temperature (6-7% per degree C) – this has led many people to postulate that atmospheric rivers in a future warmer climate will carry more water and thus produce more intense rainfall.
Continuing with our equations#
In class, we talked about the mixing ratio being approximately equal to the specific humidity and being approximately equal to 0.622(ev/p) . Lab 2-3 shows that this is true. Following Gerard, 2005 (as in the printed lecture notes), we can use this equation specifically for saturated conditions to get: $\( q_{sat}(T,z) = 0.622(\frac{e_{sat}(T)}{p(z)}) \)$
If we combine this with the ideal gas law and the hydrostatic equation, and assume that atmospheric temperature varies linearly with height, as $\( T = T_0 +\Gamma z \)\( where \)\Gamma\( is the environmental temperature lapse rate, and integrating from a height \)z\( to the surface, we get \)\( \rho q_{sat}(z) = \rho_0 q_{0 sat}\exp{(\frac{-z}{H_m})} \)\( where \)\( H_m = - \frac{b}{a \Gamma} \)\( and a = 17.67, and b=243.5\)^\circ$C, as above.
We can define condensation of water vapor as the change in saturated moisture content with time. $\( C = - \frac{d(\rho q_{sat})}{dt} ~= - \frac{\partial(\rho q_{sat})}{\partial z} \frac{dz}{dt} = - w \frac{\partial(\rho q_{sat})}{\partial z}\)$
What this says is that condensation is a function of vertical velocity, w, multiplied by the rate of change in saturated vapor pressure with height. The saturated vapor pressure decreases with height, as temperatures cool, which causes condensation as air rises.
Next, for orographic precipitation, we presume that the mountain slope is the primary source of lifting, i.e., vertical velocity. Combining the equations above, defining the mountain slope at \(\nabla z_s\), and integrating over a total vertical column of the atmosphere, we get the following:
where \(S\) is the column-integrated condensation rate (formation of water droplets)
\(\rho_0\) is inital air density (kg/m\(^3\))
\(q_{0 sat}\) is the initial specific humidity (g water/kg air)
\(\vec{u}\) is the average incoming wind speed (m/s)
\(\nabla z_s\) is the topographic slope, if we define x as the direction of the wind, dz/dx (-)
\(z_s\) is the surface elevation (must be in km if the lapse rate is in km)
\(H_m\) is the scale height for moisture as defined above (\(-\frac{b}{a \Gamma}\) (km)
and
\(\Gamma\) is the environmental lapse rate (can presume -6 \(^{\circ}\)C /km)
# Read in a .xlsx file of the height transect along a mountain range.
data = pd.read_csv('elevation-profile_Olympics.csv')
# You may also substitude with 'elevation-profile_MtRainier.csv' if you want something else.
data.head()
| lat | long | elevation_feet | distance_mi | |
|---|---|---|---|---|
| 0 | 47.70329 | -124.40506 | 130.75 | 0.00 |
| 1 | 47.70516 | -124.39059 | 161.51 | 0.69 |
| 2 | 47.70702 | -124.37611 | 175.66 | 1.37 |
| 3 | 47.70889 | -124.36164 | 210.99 | 2.06 |
| 4 | 47.71075 | -124.34716 | 240.02 | 2.74 |
These data are from the USGS Streamstats page, using their “Elevation Profile Tool” under “Exploration Tools” – after drawing a line in line with where I imagine the primary wind from an atmospheric river will hit, I scrolled down in the output and exported a .csv – (Note that I had to change their headers to not have parentheses to work with this notebook.)
We will presume that the wind is directly in line with the drawn transect, so that wind dot the gradient of elevation is the same as wind speed times the vertical change with distance here, dz/dx.
Convert to metric units, and then plot the transect data.
1 foot = 0.3048 m
1 mile = 1609.34 m
elev_m = 0.3048*data.elevation_feet
distance_m = 1609.34*data.distance_mi
fig, ax = plt.subplots(figsize=(10,4))
ax.plot(distance_m,elev_m,c='k',label='topographic profile');
ax.set_title('Elevation vs distance along the Olympic Coast', fontsize=15)
ax.set_xlabel('distance (m)', fontsize=12)
ax.set_ylabel('elevation (m)', fontsize=12);
plt.legend(loc=(0.2,-0.36), fontsize=12);
Now, we can calculate the slope along this line, and presuming we know the wind speed aligned with it, we can calculate the condensentation rate all along the line. If raindrops that condense fall out immediately and right below where they condense, we can estimate rain rates. If we do this, what physical processes are we neglecting, and how do you think we could improve on this?
# First, let's calculate dz/dx along our slope.
dx = np.diff(distance_m)
dz = np.diff(elev_m)
dzdx = dz / dx
# Note that this only gives us one value between every pair of x and z values in the plot above. For corresponding
# distance and elevation, let's calculate all the midpoints
L=len(elev_m)
x_mid = 0.5*(distance_m.values[:-1]+distance_m.values[1:])
z_mid = 0.5*(elev_m.values[:-1]+elev_m.values[1:])
#x_mid = 0.5*(distance_m[0:L-2]+distance_m[1:L-1])
#z_mid = 0.5*(elev_m[0:L-2]+elev_m[1:L-1])
# check sizes
print(len(dzdx))
print(len(x_mid))
print(L)
# Now, lets assume some typical values
ro0 = 1.2 # kg/m^3
u = 10 #m/s
qsat0 = 8 # g/kg
Gamma = -6 # degrees/km
a=17.67
b=243.5
Hm = -b/(a*Gamma)
row = 1000 #kg/m^3. density of water
# Now, combine these with the equation above and plot S/row along the line.
50
50
51
fig, ax = plt.subplots(figsize=(10,4))
ax.plot(x_mid,dzdx,c='r',label='slope');
ax.set_title('Slope vs distance from the Olympic Coast', fontsize=15)
ax.set_xlabel('distance (m)', fontsize=12)
ax.set_ylabel('slope (-)', fontsize=12);
plt.legend(loc=(0.2,-0.36), fontsize=12);
And with this, we can calculate the instantaneous condensation rates along the track#
S = ro0*qsat0*u*dzdx*np.exp(-z_mid/1000/Hm)
# Think about why it is important to divide by 1000 in the exponent here.
print(S)
[ 0.79497486 0.36997253 0.90738293 0.75331997 5.01222993
1.46712913 -3.18219638 2.73139873 12.17336383 4.63600451
-3.19858342 -2.44180397 -4.42762035 10.42584798 9.20735109
-17.77653834 -3.85678127 -0.28502599 12.01973021 12.78691437
7.27308185 3.25591325 4.84627884 1.62361168 -10.7084716
4.361424 4.78912686 -5.52117118 9.77123423 4.27913674
1.57199867 -29.40046832 18.5350356 -19.98699866 0.70464057
2.68815948 0.91297306 -17.23267986 45.62006545 4.91820355
-13.97023206 15.6963278 25.45935478 -10.28773356 1.26234412
8.2706899 -3.95374594 16.17481546 10.10257201 -14.62495011]
fig, ax = plt.subplots(figsize=(10,4))
ax.plot(x_mid,S,c='r',label='condensation');
ax.set_title('Droplet condensation vs distance from the Olympic Coast', fontsize=15)
ax.set_xlabel('distance (m)', fontsize=12)
ax.set_ylabel('condensation (g/(m^2 s)', fontsize=12);
plt.legend(loc=(0.2,-0.36), fontsize=12);
# to convert to units of rain rate in mm/hr (which make more sense)
# we divide by the density of water and convert units
rho_w = 1000 #kg/m^3
# We know that we have 1000 mm/ m and 1 kg/1000 g, and 3600 s/hr, so we end up with units of mm/hr
Rainrate = S/rho_w * 3600
# We also know that we can't get negative rain (there may be some evaporation as the parcels
# descend local places where the slope goes down, but we will neglect that for now.
Rainrate[Rainrate < 0] = 0
# the above finds all values less than 0 in the array and sets them at 0
fig, ax = plt.subplots(figsize=(10,4))
ax.plot(x_mid,Rainrate,c='r',label='condensation');
ax.set_title('Instantaneous rainrate vs distance from the Olympic Coast', fontsize=15)
ax.set_xlabel('distance (m)', fontsize=12)
ax.set_ylabel('rain rate (mm/(hr)', fontsize=12);
plt.legend(loc=(0.2,-0.36), fontsize=12);
Does this plot look reasonable to you?#
One important engineering skill is to be able to look a computer output and decide if the output is “in the ballpark” or if something went quite wrong.
If I look up tables of “probable maximum precipitation” for “mountains of western washington,” the highest number I can find is about 2 inches of precipitation in 2 hours. That corresponds to 25.4 mm/hr. Clearly, measured precipitation on the ground does not match what I just plotted.
What reasons can you think of why?#
This paper by Smith and Barstad (2004) explains the main problems with a simple linear model for orographic precipitation and develops mathematical tools to overcome them.
One simple way to improve our plot is to realize that the atmosphere responds to a smoother version of the mountain slope. If I smooth the slope over each 10 km in the horizontal, I get a smoother slope and more uniform rain rates.
# We can define a function to smooth the slope and apply it
def moving_avg(x, n):
cumsum = np.cumsum(np.insert(x, 0, 0))
return (cumsum[n:] - cumsum[:-n]) / float(n)
n=15
smootherslope = moving_avg(dzdx,n)
smootherz = moving_avg(z_mid,n)
smootherx = moving_avg(x_mid,n)
dzdx_smooth = np.diff(smootherz)/np.diff(smootherx)
x_midsmooth = 0.5*(smootherx[:-1]+smootherx[1:])
z_midsmooth = 0.5*(smootherz[:-1]+smootherz[1:])
fig, ax = plt.subplots(figsize=(10,4))
ax.plot(smootherx,smootherslope,c='r',label='slope');
ax.plot(x_midsmooth,dzdx_smooth,c='b',label='slope2');
ax.set_title('Slope vs distance from the Olympic Coast', fontsize=15)
ax.set_xlabel('distance (m)', fontsize=12)
ax.set_ylabel('slope (-)', fontsize=12);
plt.legend(loc=(0.2,-0.36), fontsize=12);
S_smooth = ro0*qsat0*u*dzdx_smooth*np.exp(-z_midsmooth/1000/Hm)
# Think about why it is important to divide by 1000 in the exponent here.
print(S_smooth)
[ 1.74911248 0.94253417 0.76007828 1.10526672 1.78702757 2.30018564
2.71623101 2.99742448 2.74277601 1.85187282 1.54998645 2.04393424
2.19701969 2.14351094 2.03453611 2.4716458 2.12250042 1.82867845
1.44433254 0.10530506 -0.40462774 -0.63393767 -1.38916525 -0.61525172
1.39260882 1.25589927 1.00209402 2.71791496 3.15785526 2.23081675
2.41623244 3.38241144 4.18818718 5.30486939 5.38407022]
Rainrate_smooth = S_smooth/rho_w * 3600
Rainrate_smooth[Rainrate_smooth < 0] = 0
# the above finds all values less than 0 in the array and sets them at 0
fig, ax = plt.subplots(figsize=(10,4))
ax.plot(x_midsmooth,Rainrate_smooth,c='r',label='rainrate');
ax.set_title('Instantaneous rainrate vs distance from the Olympic Coast', fontsize=15)
ax.set_xlabel('distance (m)', fontsize=12)
ax.set_ylabel('rain rate (mm/(hr)', fontsize=12);
plt.legend(loc=(0.2,-0.36), fontsize=12);
fig, ax = plt.subplots(figsize=(10,4))
ax.plot(distance_m,elev_m,c='k',label='topographic profile');
ax.plot(smootherx,smootherz,c='b',label='smoothed topography');
ax.plot(x_midsmooth,Rainrate_smooth*100,c='r',label='scaled rainrate');
ax.set_title('Elevation vs distance along the Olympic Coast', fontsize=15)
ax.set_xlabel('distance (m)', fontsize=12)
ax.set_ylabel('elevation (m)', fontsize=12);
plt.legend(loc=(0.2,-0.36), fontsize=12);
This still doesn’t look quite right. (At least when I go hiking in the Olympics on a rainy day, I do not find a special dry spot at exactly 30 km from the shore where no rain falls…) What else are we missing?
First, droplets do not fall from where they condense in the cloud to the ground instantaneously. It takes some time for droplets to grow and to fall. What factors do you think control this time scale? How could we add this process to our model?
Next, we have created a model that only varies in the x and z directions. What processes in the real world might complicate what is actually going on here? When might this 2-D assumption be an okay approximation of the real world? When would 3-D effects be more important?