Mountain Rainfall

Minor shifts in atmospheric circulation can have major impacts on rainfall patterns in mountainous regions. This has important implications for crops such as coffee which are grown at altitude.

Rainfall derives from uplift of moist air. In mountainous areas, uplift is forced by terrain and so depends on the horizontal wind vector \mathbf U. The dependence of rainfall on wind direction is illustrated below for a fictitious range of hills (dashed contours indicate the hills).

The “upslope” model for rainfall (P) at location \mathbf{r} is:

(1)   \begin{equation*} P(\mathbf{r}) = P_o + \rho_s  \mathbf{U} \cdot  \vec \nabla h\left(\mathbf{r}\right) \end{equation*}

\vec\nabla h is the gradient of the terrain height and \rho_s is related to the saturated water vapour density. P_0 is the background rainfall rate.

In reality, Equation (2) is a poor approximation.  Rain spreads downwind because there is a delay before raindrops can form, and another delay before fallout. Airflow dynamics also spreads the influence of an obstacle upwind. These effects can be included in a Fourier transformed version of equation (1):

(2)   \begin{equation*} P(\mathbf{r}) = P_0 + \rho_s \sum_\mathbf{k} g(\mathbf{k}) exp\left( i \mathbf{k}\cdot \mathbf{r}\right) i \mathbf{U} \cdot \mathbf{k} \tilde h(\mathbf{k}) \end{equation*}

\tilde h(\mathbf{k}) is the Fourier transform of the terrain height and g(\mathbf{k}) is a form factor which includes the effects of advection and mountain wave dynamics. Equation (2) reduces to the upslope model when g(\mathbf{k})=1.

R: Fast Fourier Transform & Rcpp

The figure illustrates orographic rainfall for a fictitious range of gaussian hills. It is very quick to compute (2) using Discrete FFT. If you use R, things can be speeded up by constructing the integrand in C++ using Rcpp before passing it to fft. This is very fast even for realistic landscapes.

C++/Rcpp implementation

R. Smith & I. Barstad, A linear theory of orographic precipitation, J. Atmos. Sci., 61, 1377–1391, (2004)




  1. Thanks for the post, really interesting… I would like to try implementing it in a landscape evolution model that I am developing.

    I was just wondering if you could send me a bit more details on how to set up the fft in R based on your Rcpp implementation. For example what are the definitions of the following parameters in the cpp file:
    – L,
    – H_a,
    – ntau1, and
    – ntau2



  2. Hi Tristan,

    L, H_a, ntau1 and ntau2 are dimensionless parameters describing the form factor g(k) of Smith & Barstad’s model.

    In terms of a (the lattice spacing of the discretised landscape), L= scorer*a, H_a = vertical scale in units of a, ntau = U/a*tau where tau are timescales defined in Smith and Barstad. U is the wind speed.

    the mountain wave dynamics is controlled by the scorer parameter (Brunt–Väisälä frequency of the stable atmosphere divided by U).

    The R implementation consists of

    (1) fft on the gradient of the landscape matrix S(k,l)
    (2) multiply by form factor transMatC_Re + i*transMatC_Im
    (3) inverse fft


Leave a Comment

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.