Evaluation of a High Resolution Wind Model Over 

a Complex Terrain Surface 

by Sam S. Chang, Dennis M. Garvey, 

Chatt C. Williamson, and Giap Huynh 

ARL-TR-4773 April 2009 

Approved for public release; distribution unlimited. 
NOTICES 

Disclaimers 

The findings in this report are not to be construed as an official Department of the Army position 

unless so designated by other authorized documents. 

Citation of manufacturer’s or trade names does not constitute an official endorsement or 

approval of the use thereof. 

Destroy this report when it is no longer needed. Do not return it to the originator. 
Army Research Laboratory 

Adelphi, MD 20783-1197 

ARL-TR-4773 April 2009 

Evaluation of a High Resolution Wind Model Over 

a Complex Terrain Surface 

Sam S. Chang, Dennis M. Garvey 

Chatt C. Williamson, and Giap Huynh 

Computational and Information Sciences Directorate, ARL 

Approved for public release; distribution unlimited. 
ii

REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the 

data needed, and completing and reviewing the collection information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing the 

burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 

Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid 

OMB control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY)

April 2009 

2. REPORT TYPE 

Final 

3. DATES COVERED (From - To) 

5a. CONTRACT NUMBER 

5b. GRANT NUMBER 

4. TITLE AND SUBTITLE 

Evaluation of a High Resolution Wind Model Over a Complex Terrain Surface 

5c. PROGRAM ELEMENT NUMBER 

5d. PROJECT NUMBER 

5e. TASK NUMBER 

6. AUTHOR(S) 

Sam S. Chang, Dennis M. Garvey, Chatt C. Williamson, and Giap Huynh 

5f. WORK UNIT NUMBER 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

U.S. Army Research Laboratory 

ATTN: AMSRD-ARL-CI-ED 

2800 Powder Mill Road 

Adelphi, MD 20783-1197 

8. PERFORMING ORGANIZATION 

 REPORT NUMBER 

ARL-TR-4773 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S ACRONYM(S) 

11. SPONSOR/MONITOR'S REPORT 

 NUMBER(S) 

12. DISTRIBUTION/AVAILABILITY STATEMENT 

Approved for public release; distribution unlimited. 

13. SUPPLEMENTARY NOTES 

14. ABSTRACT 

A complete and comprehensive description of the high resolution wind (HRW) model is presented. The HRW model has been 

in development at the U.S. Army Research Laboratory, formerly U.S. Army Atmospheric Sciences Laboratory, since 1978. 

This model is a two-dimensional, diagnostic atmospheric surface-layer wind model with a horizontal grid spacing of the order 

of 100 m over a domain of about 5 by 5 km. It uses Gauss’ principle of least constraint and a direct variational relaxation 

method to adjust an initially uniform wind field to conform with topography, mass conservation, and buoyancy forces. A 

distinctive feature of the model is the use of a non-orthogonal, terrain following, warped coordinate system. A valuable 

observational dataset of surface wind has been provided from the field study of Meteorology and Diffusion Over Non-Uniform 

Areas (MADONA) at Porton Down, Salisbury, England during September and October 1992. Using the MADONA data, a 

critical evaluation for the HRW model for 39 cases has been carried out. The results of this study are presented, indicating both 

the range of validity and the limitations of the HRW model. 

15. SUBJECT TERMS

High resolution wind (HRW) model, atmospheric surface layer 

16. SECURITY CLASSIFICATION OF: 

19a. NAME OF RESPONSIBLE PERSON

Sam S. Chang 

a. REPORT 

Unclassified 

b. ABSTRACT 

Unclassified

c. THIS PAGE 

Unclassified

17. LIMITATION 

 OF 

 ABSTRACT 

UU 

18. NUMBER 

 OF 

 PAGES 

45 19b. TELEPHONE NUMBER (Include area code) 

(301) 394-1988 

Standard Form 298 (Rev. 8/98) 

Prescribed by ANSI Std. Z39.18 
iii

Contents 

List of Figures iv

List of Tables iv

Acknowledgments v

1. Introduction 1

2. Description of the HRW Model 2

2.1 Gauss’s Principle of Least Constraint ............................................................................ 2 

2.2 Warped Coordinate System............................................................................................ 4 

2.3 Empirical Surface Layer Profiles ................................................................................... 9 

2.4 Computational Aspects and Algorithm ........................................................................ 11 

2.4.1 Data Input and Model Setup............................................................................. 12 

2.4.2 Relaxation Scheme ........................................................................................... 14 

2.4.3 Data Output ...................................................................................................... 15 

3. Results of the Model Evaluation 16

3.1 The MADONA Data .................................................................................................... 16 

3.2 The Model Simulations ................................................................................................ 17 

3.3 The Results of Linear Regression Analysis ................................................................. 19 

3.4 Discussion .................................................................................................................... 28 

4. Conclusions 30

5. References 32

Distribution List 34 
iv

List of Figures 

Figure 1. The terrain following, warped and non-orthogonal coordinate system with three base 

vectors 12 3 ( , , and ) aa a

  which is defined by equations 11, 12, and 13, respectively. V



 is a 

wind vector. A computational volume element with top (HMNP) and bottom (OEFG, 

shaded) is also shown, see text for details. ............................................................................... 5 

Figure 2. Computational flow chart for the HRW model. This chart illustrates the three essential 

stages and their major routines. .............................................................................................. 12 

Figure 3. Locations of MADONA measurements sites (Ml through M15) plotted on terrain 

contours................................................................................................................................... 17 

Figure 4a. Scatter diagrams for simulated wind direction (y, degree) versus observed values (x, 

degree). The straight lines result from the linear regression analysis, equation 42................ 20 

Figure 4b. Scatter diagrams for simulated wind speed (y, ms

–1) versus observed values (x, ms–1). 

The straight lines result from the linear regression analysis, equation 42.............................. 21 

Figure 5. Same as figure 4a except for three wind direction categories as defined by equation 

43a, b, c................................................................................................................................... 24 

Figure 6. Same as figure 4b except for three-wind speed categories as defined, by equation 43d, 

e, f. .......................................................................................................................................... 25 

Figure 7. Same as figure 4 except for three stability categories as defined by equation 44a, 

b, c........................................................................................................................................... 27 

List of Tables 

Table 1. A list of input and output data for the HRW model....................................................... 13 

Table 2. The 39 cases for the HRW simulations. Q is defined by equation 24a in the text...... 18 

Table 3. The value of n (data pair number), a (slope), b (intercept), r (correlation coefficient), 

and S (standard deviation) from the linear regression equation 42, for wind direction and 

wind speed, respectively. ........................................................................................................ 22 

Table 4. Same as table 3 except for three wind direction categories as defined by equation 43 

(a, b, c) and three wind speed categories as defined by equation 43 (d, e, f). ........................ 23 

Table 5. Same as table 3 except for three atmospheric stability categories defined by equation 

44a, b, c................................................................................................................................... 28 
v

Acknowledgments 

We appreciate Mr. Ronald Cionco, our now retired colleague, who played a leading role in the 

early development of the model, and continued development and application of the model for 

various scenarios. Dr. Jon Mercurio, also retired colleague, contributed to the construction of the 

model and this manuscript. Dr. James Cogan and Dr. Douglas Brown offered administrative 

support for this technical report. 

Miss Patricia Cordes, a student at the Academy of the Holy Cross, who participated in the 

Science and Engineering Apprentice Program sponsored by the U.S. Army Research Laboratory 

in Summer 1999, contributed much of the linear regression analysis and graphing. Mr. John 

Byers helped check the data files used in this study. Mr. Ron Meyers provided useful 

discussions on the model. Mr. Richard Hasenauer helped with the figures. Ms. Deborah 

Lehtinen has typed the manuscript with great skill. 
vi

INTENTIONALLY LEFT BLANK
1

1. Introduction 

During the last three decades meteorological researchers at the U.S. Army Research Laboratory 

(ARL), formerly the Atmospheric Sciences Laboratory (ASL), have conceived, developed, and 

used a high resolution wind (HRW) model. The HRW model is a micro-alpha scale (Orlanski, 

1975), two-dimensional, diagnostic model that simulates wind and temperature fields in the 

atmospheric surface layer, taking into account both complex terrain topography and thermal 

structure over a limited area. The computational domain size can range from a 2 km square to a 

20 km square with grid resolutions varying from 40 m to 400 m, respectively. The vertical 

thickness of the computational layer is designed to be one tenth the magnitude of the grid size. 

A typical grid size of 100 m, therefore, produces simulated fields at a level 10 m above the 

surface. 

The HRW model was originally formulated by Ball and Johnson (1978) who described its 

theoretical basis and computational structure in a technical report submitted to ASL. This 

geographically re-locatable model was designed to be incorporated in the U.S. Army 

Experimental Prototype Automatic Meteorological System for the estimation of the surface layer 

wind field at sub-mesoscale resolution over a limited area in broken topography. Later the HRW 

model was further developed as a stand-alone model, as well as an integral component of a 

hierarchy of nested meso- and micro-meteorological models (Cionco, 1987). A distinctive 

feature of the HRW model is that it has adopted a warped coordinate system to address the 

intimate interaction between the surface layer and the variable ground features. As far as the 

authors are aware, this approach is unique; virtually no meso- or micro-scale models have been 

formulated with this approach. The HRW model has been tested and used for a variety of 

applications, e.g., Weber et al. (1995), Thykier-Nielsen et al. (1995), Cionco and Byers (1995, 

1997), Cionco (1998), Cionco et al. (1998), Cionco and Ellefsen (1998). Although the model has 

been tested and used for a long time, it has not been thoroughly evaluated due to the fact that 

adequate observational data at a resolution of 100 m or so are extremely scarce. Fortunately, the 

multinational, high-resolution field study of Meteorology and Diffusion over Non-Uniform 

Areas (MADONA) during September and October 1992 has provided a valuable observational 

dataset. Cionco et al. (1999) have given a comprehensive description of the MADONA project, 

including the dataset. The MADONA field study was designed and conducted for, among other 

purposes, high-resolution meteorological data collection in an effort to obtain terrain-influenced 

meteorological data. Thirty-one days of meteorological data were collected. High-resolution 

and standard micro-scale, boundary layer, and synoptic meteorological sensors including 15 

wind speed/direction sets were deployed over the MADONA topography (a 9 km by 7.5 km 

area). This well-documented database is suitable and valuable for the evaluation and validation 

of the HRW model. 
2

The objectives of this report are twofold. The first objective is to provide a concise description 

of the HRW model as systematically and completely as possible. The lengthy paper by Ball and 

Johnson (1978) was not published in the open literature and is not readily available for interested 

readers. A limited number of papers presented in a book (Cionco, 1985) or at conferences (cited 

above) have separately described only gross features of the model without providing necessary 

details. Section 2 provides the basic theory, warped coordinate system, mathematical 

formulation, and the computational algorithms of the model. The second objective is to present 

the results of a critical evaluation of the model using the MADONA database (section 3). 

Finally, the last section (section 4) gives conclusions from the present study of the HRW 

model, and relates these results to a companion evaluation of the model in a report by 

Williamson et al. (2005). 

2. Description of the HRW Model 

2.1 Gauss’s Principle of Least Constraint 

The HRW model uses Gauss’s principle of least constraint directly rather than the more 

commonly employed Newton’s laws of motion. A discussion of this principle applied to systems 

of point particles can be found in Lanczos (1966). As applied to non-viscous, incompressible 

fluids, the principle can be expressed as 

  1 2

0 , 2    A g d pA d                (1) 

where  denotes a variation of the two integrals in square brackets. The vector A

 is the 

acceleration, dV dt / ,  where V



 is the velocity vector; and g



 is the acceleration of gravity, the 

only external force considered. The two integrals above are over a material volume () and its 

boundary surface (), respectively. The symbols  and p are the fluid density and pressure 

respectively. For the atmosphere,  and p are related by the equation of state for an ideal gas, 

p  RT, (2) 

where R is the specific gas constant for air, and T is the absolute temperature. The more 

conservative potential temperature  is related to T as 

  /

/ ,

R cp

  Tp p ref (3) 

where pref is the reference pressure taken as 100 kPa and cp is the specific heat at constant 

pressure. 

Equating the variation of the expression in brackets to zero implies that the air motion takes 

place in such a way as to minimize constraint forces arising solely from kinematic conditions. In 
3

this case the constraints are the boundary conditions and conservation of mass, which for an 

incompressible fluid is expressed as 

  V 0.  (4) 

This constraint results in the pressure gradient acting as a force of constraint necessary to enforce 

incompressibility. Ball and Johnson (1978) discuss this result in greater detail. 

In order to account for atmospheric thermal effects, additional assumptions are needed. First, the 

variables in the equation of state are decomposed as 

0 000              , ,,, p p pT T T (5) 

where the variables with zero subscripts denote ambient or mean values while the same variables 

with primes represent departures from the ambient. Secondly, using the Boussinesq 

approximation, equation 1 can be approximated by 

 2

0

1 0

2    A bg d p A d                 (6) 

b     / /,   0 0    (6a) 

where , b g



 the effective external acceleration, is the buoyancy acceleration, and b is the 

buoyancy parameter defined by equation 6a. A detailed discussion of the Boussinesq 

approximation is provided by, e.g., Stull (1988, Chapter 3). Likewise the effective pressure in 

equation 6 is the fluctuating pressure ( ) p , the departure from the ambient pressure (p0). It is 

assumed that both p and p0 satisfy the hydrostatic equation 

0 0 dp g dz dp g dz     , ,   (7) 

where z denotes vertical height. Thirdly, similar to equation 6 a variational expression can be 

applied to the potential temperature field 

2

1 0, 2

d d

dt

  

        

       (8a) 

where the adiabatic condition has been assumed. Equation 8a implies that 

0. d V

dt t

       



 (8b) 

For a steady state, equation 8a is equivalent to 

  

2

  V d 0.           (8c) 
4

This equation is needed in order to make the temperature field and the buoyancy forces 

consistent with the wind field. 

For the application of the basic variational principles, equations 6 and 8c, to a single surface 

layer, further simplification is required. This single surface layer is of constant thickness in a 

direction normal to the terrain surface and follows the warped ground surface. The second 

(surface) integral in equation 6 contains the pressure fluctuation p . This p term is extremely 

difficult to estimate. A practical alternative is to reformulate equation 6. With the help of 

Gauss’s theorem, equation 6 can be rewritten for the single surface layer as 

   2

0

1 0, 2 A A bg bg p A A d sn s n sn                          

where the subscripts s and n denote the surface parallel and surface normal component, 

respectively. The integrand in the above equation can be rewritten as 

    2 2

0

1 2

2

 A A bg bg bA g p A p A s n sn n n ns

              

       

One part of the second-to-last term in the above expression, ( )  pAn

  can be used to cancel the 

last term within the brackets, 0 ( ) n n  bA g

 

 since / n p   n g  as implied by equation 7. The 

absence of the cross product of A n

 and n b g suggests that n b g may be ignored. The other part 

of the second-to-last term, ( ) n p A  

 as well as the last term, ( )s  p A

  can be neglected because 

p is usually small Consequently equation 6 is further simplified as 

 

2

0 A bg d s              (9) 

where only the surface parallel component of buoyancy s b g

 remains, and the multiplicative 

constant (0/2) has been neglected. The surface parallel simplification in equation 9 assumes that 

for the surface layer the acceleration normal to the terrain surface is not affected by the normal 

buoyancy force ( ) n bg . 

2.2 Warped Coordinate System 

A terrain-following, warped, non-orthogonal coordinate system is employed in the HRW model 

in order to account for a complex and varying terrain surface as accurately as possible. Figure 1 

is a diagram of the warped coordinate system. The upper portion of figure 1 indicates a 

Cartesian coordinate system with respect to a horizontal reference plane in which , , i j

  and 

k



are three unit vectors along the three orthogonal directions (x, y, and z), respectively. Any 

position on the terrain surface is given by the vector 

 (, ) , r xi y j h x y k       (10) 
5

where h(x, y) is the surface elevation. 

Figure 1. The terrain following, warped and non-orthogonal coordinate 

system with three base vectors 12 3 ( , , and ) aa a

   which is defined 

by equations 11, 12, and 13, respectively. V



 is a wind vector. 

A computational volume element with top (HMNP) and bottom 

(OEFG, shaded) is also shown, see text for details. 

The lower portion of figure l shows three base vectors, which are defined by 

2

1 11 ,, 1 x x x

r h a i hk h a a h

x x

 

     

 

     (11) 

2

2 22 ,, 1 yy y

r h a j hk h a a h

y y

 

     

 

     (12) 

  2 2

3 3

1 , 1 , 1. x y xy a n hi h j k a h h a

a

             (13) 

These three vectors 12 3 aa a , , and   establish a warped coordinate system. Notice that 1 a

 and 2 a



are parallel to the surface, but are not unit vectors. The vectors 1 a

 and 2 a are not perpendicular 

to each other, and hence, the warped coordinate system is not an orthogonal coordinate system. 
6

The vector 3 a

 is a unit vector and is perpendicular to the terrain surface 1 2 ( and ) a a  . The volume 

of the parallelepiped formed by these three base vectors is 

aaa a 123    .    (14) 

The angle between 3 a

 and the vertical z axis is  as indicated in figure l, 

cos 1/ .   a (15) 

It is useful to define the system of reciprocal base vectors denoted by a superscript 

    1 2

2 3

1 / 1 y xy x a a a a h i hh j h k

a

        

     (16) 

    2 2

3 1

1 / 1 xy x y a a a a hh i h j h k

a

         

     (17) 

  3

12 3

1 / . x y a a a a hi h j k a

a

        

      (18) 

The reciprocal vector 3

a

 is identical with 3 a or n

, the unit vector normal to the terrain surface. 

The reciprocal vectors 1

a

 and 2

a

 are parallel to the surface, but distinct from 1 a

 and 2 a . Also 

1

a

 and 2

a

 are not perpendicular to each other. The system of reciprocal vectors facilitates 

operations because 

, 1,2,3, s

r rs a a rs      (19) 

where rs = 1 if s = r, and rs = 0 otherwise. 

A three-dimensional wind velocity vector in the atmospheric surface layer can be expressed as 

12 3

123 . i V Va V a Va Va    i

   (20) 

where 1 2

1 2 V a OA V a OC   ,     , and 3 V a OD 3 

  as illustrated in figure 1. The vector OB

 is the 

component of the wind velocity along the terrain surface. The sum convention for repeated index 

(i) is used in equation 20 and hereafter. Notice that 2 12 22 32 () ( ) ( ) ( ) VV V V  due to nonorthogonality of the coordinates. The component V 3 is called the surface normal velocity 

component or the impaction vertical motion. Note that V 3

 is not the vertical velocity (w) in 

terms of a Cartesian coordinate system. It is not difficult to prove that 

  1 2 1 1 V h u hhv hw y xy x a

     (21a) 

  2 2 1 1 V hhu h v hw xy x y a

        (21b) 
7

3 1 V hu hv w x y

a       (21c) 

where u and v are the wind velocity components along x and y direction, respectively (figure 1). 

The surface kinematic boundary condition is usually expressed by 

 ( ) (, ) wz h V hxy     (22a) 

with the help of equation 21c, the above condition leads to 

V 3

(z = h) = 0 (22b) 

That means the surface kinematic condition imposes zero impaction vertical motion, but not the 

Cartesian vertical velocity (w) at the terrain surface. As shown in figure 1, the acceleration of 

gravity is 

        g gk g n S g g cos sin     n s        (23) 

where n g

 and s g

 denote the surface normal and surface parallel component of g



, 

respectively. The term S



 is a unit vector along the OB



 direction in figure 1. As discussed 

earlier, s bg is used for the calculation of buoyancy acceleration in equation 9. Besides the 

buoyancy parameter (b), a surface heating (dq) is introduced in the model. From the First Law 

of thermodynamics for an ideal gas, the heating can be written as 

dq = cp dT – dp/ = cp T(d/) 

where the definition of , equation 3, has been used, see, e.g., Wallace and Hobbs (1977). 

Hence for the surface layer, dq can be related to the buoyancy parameter (b) as defined by 

equation 6a 

0

0

, sfc

p sfc p sfc dq c T c T b

 



        

 (24) 

where the subscript “sfc” denotes the surface layer. In order to evaluate the integrals of 

equations 8c and 9, we use a computational volume element for the single surface layer as 

sketched in figure 1. This volume element, also called local flux box, consists of the layer 

between the terrain surface, ( )  and a constant normal height (d) with lateral lengths l1 and l 2

as indicated in figure 1 

2

1 1 1  OE GF la l h    x  (25a) 

2

2 2 1 .    OG EF la l hy  (25b) 
8

where l OE E F    is the side length of the square which is the projection of the volume 

element on the horizontal reference plane. The volume element is not rectangular and its 

volume () is 

   2

1 11 2 2 2 3      l a a l a a da l d a // .    

   (26) 

For reference, the lateral faces of the local flux box are numbered counterclockwise (1), (2), 

(3), and (4) as shown by EFNM, FGPN, GOHP, and OEMH, and (5) and (6) denote top 

(HMNP) and bottom (OEFG) faces, respectively in figure 1. The vector expressions of the 

six surface areas are 

  1

1 2 22 3 1 1

   S l a a da ld aa   /      

    (27a) 

  2

2 3 1 11 2 2

    S d a l a a ld aa   /      

    (27b) 

  1

3 2 22 3 3 3

    S l a a da ld aa   /      

    (27c) 

  2

4 3 1 11 4 4

    S da l a a ld aa   /      

    (27d) 

   2 3

5 1 11 2 2 2 5 5

   S l a a l a a l aa   / /      

    (27e) 

   2 3

6 1 11 2 2 2 6 6

    S l a a l a a l aa   // ,      

    (27f) 

where 1 2

a a, ,   and 3 a

 have been defined by equations 16, 17, and 18 respectively. The 

subscripts: (1 through 6) in equation 27 mean that all terrain dependent quantities must be 

evaluated using representative (average) values on that surface. The surface area vector 

  S k k ( 1, 6)   is in the outward normal direction of each face. Using the divergence theorem 

for the volume element, the condition of mass conservation equation 4 can be approximated by 

    12 3

123

1 1 0.

k k k k

V V S Va Va Va S

 

                       (28a) 

From equations 19 and 27 the above condition yields 

            3 3 12 12

5 6 1 23 4

, d d aV aV aV aV aV aV W

l l

          (28b) 

where W is only a notation to be used later. The term 3

6 ( ) aV in equation 28b vanishes due to 

equation 22b. The surface normal velocity (V 3

) in the surface layer is of order (d/l) < 0.1 with 

respect to surface parallel components (V 1

 and V 2 term) as implied by equation 28. Using 
9

equations 19, 20, and 27, and dyadic (tensor) analysis (Tai, 1992, for example) the acceleration 

vector for the local flux box can be written as 

    

  

1 11 21 31

123

1

12 22 32

123 2

11 21 31

123 3

12 22 32

123 4

13 23

1

ij ij 1

ij ij

k k

A VV V V a a V V a a S

l a V V aa V V aa V V aa

V V aa V V aa V V aa

V V aa V V aa V V aa

V V aa V V aa V V aa

l V V aa V V a

d





     

  

     

     

     

     

 

          









 3 3

2 3 5

13 23 33

123 6

a V V aa

l V V aa V V aa V V aa

d

    

        

 



 (29) 

With the help of equations 11, 12, and 13 A



 can also be expressed for the Cartesian coordinate 

system as 

          

           

    

   

  

11 12 11 12 1

1 23 4 5

21 22 21 22 2

1 23 4 5

11 21 12 2 2

1 2

11 21 12 2 2

3 4

1 2

5

1

1

1

,

xy x y

xy x y

x y

A i VV a VV a VV a VV a VW

l a

j VV a VV a VV a VV a VW l a

k a VVh VVh a VV h VV h

l a

a VVh VVh a VV h VV h

W Vh Vh

   

   

        

        

     

 



 (30) 

in which all of the terms of V 3

Vm (m = 1, 2, or 3) have been neglected since V 3 « V 1 or V 2 and 

equation 28b has been used. V3

 can be large only if hx and/or hy in equation 21c is large. 

2.3 Empirical Surface Layer Profiles 

Accurate calculation of the integrals in equations 8c and 9 for a volume element requires use of 

surface normal profiles of temperature and wind velocity in the surface layer. The HRW model 

assumes that the following empirical profiles in the z direction may still be applied locally in the 

n direction (figure 1). The model has adopted profiles for wind speed and potential temperature 

as follows: 
10

() ,

m

r

r

z Vz V

z

      

 (31) 

and 

(z) = r + (z – zr) , (32) 

where Vr and r denote the wind speed and the potential temperature at a reference height (zr) 

which can be the top of the volume element. From these equations note that the potential 

temperature is assumed linear with height,  being the gradient, and that the wind speed follows a 

power law (e.g., Arya, 1999). The value of the exponent m can be estimated from surface layer 

similarity theory, which defines a non-dimensional wind shear 

  . m m

z kz dV

L u dz  



        (33) 

Hence 

 

0 0

( ) ln 1 ln ln , m m

u u z z zz Vz d

kz L kz L                            (34) 

where z0 is the roughness length, L is the Obukhov length, and m is a newly defined function; 

see Garratt (1992) for details. Panofsky and Dutton (1984) have shown 

0

.

ln

m

m

z

L

m

z z

z L





             

 (35) 

The empirical wind shear and wind speed profiles of Businger et al. (1971) have been used by 

the HRW model. It is impractical, however, to use their original formulation since L and hence 

the non-dimensional height z

L      

  are usually unknown. Instead, an empirical expression for 

m in terms of the Richardson number (Ri) can be used. For example, Ball and Johnson (1978) 

used the following expression: 

 

 

 

 

1

4

1

1

1

1 12 ; 0.0

1 3 ; 0.0 0.03571

0.88 1 6 ; 0.03571 0.1246

1.75 1 4 ; 0.1246 0.25

m Ri Ri

Ri Ri

Ri Ri

Ri Ri

 







 

  

  

  

 (36) 
11

in which 

Ri = (Ri)b + (Ri)c , (37) 

where the bulk Richardson number, (Ri)b, is 

 0 

2

0

() . b

r

g z

Ri

V

 



       

(38) 

(Ri)c in equation 37 represents the Richardson number for curved flows. This additional 

Richardson number was originally suggested by Bradshaw (1969) for curvature effects 

of curved flows, 

( ) 2 (1 ) , . c

V dV Ri S S S

R dR

           (39) 

Here R is the radius of streamline curvature. If the average surface normal acceleration An in the 

layer is used and is assumed to be purely centripetal (–VV/R), and the average value of V(dV/dz)

of equation 31 is used, then S can be approximated by 

2

2 . n r

r

A z S

V  

2.4 Computational Aspects and Algorithm 

The computational aspects of the HRW model are illustrated in figure 2, a simplified 

computational flow chart. As indicated in figure 2 there are three main stages: (1) data input 

and model setup; (2) relaxation scheme through direct variation; and (3) model data processing, 

output and plotting programs. These primary stages are discussed in the following sections. 
12

Figure 2. Computational flow chart for the HRW model. This chart 

illustrates the three essential stages and their major routines. 

2.4.1 Data Input and Model Setup 

To initiate a model run, necessary input data are required to provide control and driving 

parameters. Table 1 (upper part) lists the required input data. First, there are a number of 

control parameters which include the setting of the computation domain, grid system and 

resolutions, important constants, data handling and storage, etc. Secondly, the computation 

requires input of digitized terrain data at every grid point, h(i, j) i, j = 1, 2, 3, - - , N. The input 

file of terrain height for the major routine (Read-Terrain, in figure 2) should be properly 

formatted. The additional vegetation data, i.e., the vegetation element heights, hveg(i, j), are also 

needed. In such a case, the effective terrain heights, heff (i, j), will be calculated as 

 ( , ) ( , ) 0.70 ( , ). eff veg h i j hi j h i j   

An empirical expression for the roughness length (z0) over vegetation has been used 

0 ( , ) 0.014 0.14 ( , ), veg z ij h ij   

where z0 is in meters, see Cionco (1983). 
13

Table 1. A list of input and output data for the HRW model. 

 Input Category Notation (unit) Remark 

1 Control parameters D(x, y), Computation domain

 And flags x, y resolutions typically x = y = 100 m 

 grid system, relaxation 

 processes, etc 

 2 Digitized terrain data h(i, j), (m) i,j = 0, 1, 2, ...,N (N = 50) 

 3 Vegetation data hveg (i, j), (m) Vegetation element height 

 4 Surface station observation Vsfc

 (at 10 m), (ms–1) Constant value within D

 Tsfc (at 10 m), (K) 

 p(3), T(3), h(3) pressure, temperature and 

 5 Upper air sounding p(2), T(2), h(2) geopotential height at 

 p(1), T(1), h(1) three levels, e.g., at 700, 850, mb 

 and at the surface 

 Output 

 1 Calculated wind field V



(i, j) (ms–1) defined by (21) 

 2 Calculated potential sfc(i, j) (K) calculated by (24b) 

 temperature field 

 3 Richardson number Ri(i, j) defined by (37) 

 4 Friction velocity 1 u i j ms ( , )( ) 

 calculated by (34) 

 5 Wind profile exponent m calculated by (35) 

As indicated in table 1, a minimal set of meteorological data is also required to initialize the 

model. The initial wind field ( ) Vsfc

 is obtained either by local meteorological observation of 

wind speed and direction (at 10 m height) or by an estimation of the geostrophic wind near the 

surface. The model assumes a homogeneous initial wind field, i.e., the initial velocity , Vsfc

 is 

constant throughout the computational domain. On the other hand, the initial potential 

temperature, sfc(i, j), varies within the domain. To calculate the initial field of (i, j), the 

following algorithm has been adopted. 

(a) Using the upper air sounding as the input data, the ambient potential temperature (0) as 

defined by equation 5 is derived by extrapolating the background potential temperature profile to 

the surface. For example 

   850 700 

0 850 850

850 700

, sfc h h

h h

        (40) 

where 850 and 700 are the potential temperature at the 850 and 700 mb, respectively. Further, 

hsfc, h850, and h700 refer to the geopotential height at the surface, 850 and 700 mb, respectively, 

and hsfc is also the altitude of the sounding station. 
14

(b) Also using the upper air sounding and the terrain height data, h(i, j), the surface pressure at 

each grid point is calculated by a linear extrapolation, e.g., 

  

850

850 700

850 700

pi j hi j h ( , ) 850 ( , ) , h h

   

where the p(i, j) has the unit of mb. 

(c) To calculate the surface heating (Q) as an input parameter from equation 24 

 0 

0

, sfc

sfc

p

dq Q T

c

 



   (24a) 

where the surface temperature (Tsfc) and the surface potential temperature (sfc) are obtained from 

observation of a surface meteorological station. 

(d) It is assumed that the value of Q from a single (station) observation is representative for the 

whole domain. This assumption means the right hand side of equation 24a can be applied to any 

grid point (i, j). Thus 

0.286

0

0 0

(, ) 1000 (, ) (, ) (, ) , (, ) (, ) sfc

sfc

i j ij ij Q ij Q T i j pi j  

                

 (24b) 

where the ambient background potential temperature at a gridpoint (i,j) can be extrapolated from 

the sounding data, equation 40. 

  850 700

0 850 850

850 700

(, ) (, ) . i j hi j h h h

        (24c) 

Lanicci (1985) has discussed the calculation of the  Q and has suggested some modifications. 

For a location where hsfc > h850 , we will use the sounding data at higher levels, e.g., 500 and 

700 mb instead of at lower levels, e.g., 700 and 850 mb. 

2.4.2 Relaxation Scheme 

As indicated by figure 2, the central stage of the HRW model computation is the relaxation 

scheme, which is the workhorse of the model. The scheme implements the calculation of the 

integrals, equations 8c and 9 as well as the variational technique to find the minimum of 

equation 9. As mentioned in section 2.2, the evaluation of the integral in equation 9 adopts 

computational elements called flux boxes. Consequently, the total constraint integral of 

equation 9 can be expressed as a sum over all flux boxes in the modeled area 

 

2

, , ,

. T ij s i j ij ijR R A bg      (9a) 
15

For variational purposes we consider the Ri,j, called residues, as functions of only the two 

velocity components (V1

, V2)i,j at the grid point (i, j). Denoting the two velocity components at 

each grid point by , ( 1, 2) k V k i j  , the unit vector of global steepest descent in this twodimensional velocity space is written as 

,

,

, 1

2 2

,

, ,

.

i j

k

k i j

i j

i j

k

ij k i j

R

V

n

R

V



  

                     

 

 (41) 

The steepest descent method computes a velocity correction , ( ) k Vi j proportional to the steepest 

descent unit vector for each relaxation sweep over the entire grid. This velocity correction has 

been chosen as 

, ,

k k   V C V N Mn ij r r ij  

where Cr is a selected fraction of a velocity scale, Vr. Currently Cr = 0.02 and Vr is the constant 

initial wind speed or 2 ms–1, whichever is greater. N and M are the grid dimensions, currently 

N = M = 51. All velocity corrections are applied simultaneously at the end of each sweep. 

This implies the velocity field after each relaxation step is 

11 1 22 2

,, , ,, , , . VV V VV V ij ij ij ij ij ij    

This relaxation algorithm apportions larger corrections to regions of the grid where the constraint 

integral is most sensitive to changes in the velocity field. The relaxation scheme calculates the 

total constraint RT at each relaxation sweep and saves the wind field for the minimum value of 

RT achieved. Each relaxation sweep also adjusts the potential temperature field by use of 

equation 8c. Consequently each relaxation sweep yields altered potential temperature and 

velocity fields until equation 9 is satisfied. 

2.4.3 Data Output 

The final stage of the model computation produces various computational results. As shown in 

table 1, the output data comprises at least five categories. The most important output is, of 

course, the surface layer wind field. It should be stressed that the calculated windfield refers to 

the three components in the terrain following warped coordinate system since the model is 

formulated according to such a unique coordinate system. It is, however, easy to transform these 

three components (V1

, V2, V3) into Cartesian wind components (u, v, w) as indicated by 

equation 21. The second category of output is the calculated field of surface potential 

temperature which can also be transformed into surface temperature easily. Finally the other 

three categories of output are related to turbulence characteristics (u and Ri) and the exponent 
16

(m) of the power law of wind profiles in the surface layer. During the past decade several 

versions of algorithms and programs for output displays and plotting have been tested and 

established. Those aspects related to the software developments will not be discussed 

in this report. 

3. Results of the Model Evaluation 

3.1 The MADONA Data 

The micrometeorological data used in this evaluation study are based on measurements during 

the field study of MADONA in September 1992. Cionco et al. (1999) have provided a 

comprehensive overview of the MADONA project including its database. The MADONA field 

site is located at Porton Down, near Salisbury, United Kingdom. From 14 to 23 September 1992, 

there were ten days of intensive measurements. Extensive high-resolution micrometeorological 

data were collected to complement the diffusion data. Fifteen sets of standard wind speed (three 

cup anemometer) and wind direction (vane anemometer) sensors were deployed at 14 locations 

by the UK Chemical and Biological Defence Establishment (CBDE, now part of the Defence 

Evaluation and Research Agency). These wind sensors were calibrated before and after the field 

measurements. Figure 3 illustrates the locations of these wind measurement stations (Ml through 

M15) plotted on terrain contours. All wind sensors except M9 were mounted at 10 m above the 

ground. The M9 station was installed at 30 m above the ground at the same tower as M8 to 

provide dual measurements at the same location. The area covered by figure 3 is 5 by 5 km, 

which is the domain for the HRW simulations. The wind data from Ml and M3 are excluded for 

model evaluation because their locations are very close to the model boundary, and the data from 

M8/9 are not used due to technical reasons. Altogether the measured wind data from eleven 

locations (Station M2, M4–7, and M10–15) are used for the comparison between the 

measurements and the HRW simulations. All the wind data, as well as other MADONA data 

including the area terrain data, have been officially documented on CD-ROM. The CD-ROM 

includes the averages of wind speed and direction for 5 min, 10 min, 30 min, and 60 min, among 

other statistics. The 5 min averages of wind measurements at 10 m above the ground from these 

11 locations are those used for the present study. Unfortunately, there were no corresponding 

temperature sensors at the same 11 stations provided by CBDE due to financial constraints. 

The CD-ROM is available on request from the RISO National Laboratory, P.O. Box 49, 

Roskilde, Denmark. 
17

Figure 3. Locations of MADONA measurements sites (Ml through M15) plotted 

on terrain contours. 

As shown in figure 3, the topography at the MADONA site is gently rolling terrain with a ridge 

running southwest to northeast (approximately 230°–50°) with higher terrain at each end. Total 

relief for the domain is about 100 m from 80 m northwest to 180 m near Tower Hill. In such a 

small area, however, the wind distribution inhomogeneity is significant, as shown by Cionco et 

al. (1999, their figures 7 and 8). Generally, the wind speeds over the ridge (e.g., M14) are 

stronger, while over the lower land (e.g., M13) they are weaker. 

3.2 The Model Simulations 

During the ten days of intensive measurements from 14 to 23 September 1992, a series of 

boundary layer soundings were launched for each date from the middle of the MADONA 

experimental area, as described by Cionco et al. (1999). Also, traditional upper air soundings 

from the UK Meteorological Office Larkhill station and surface meteorological observations 

from CBDE’s Met office were available during the experiment. These sounding and surface 

observation data (surface temperature at 10 m AGL in particular) provide the necessary inputs 

for the HRW simulations. The wind data measured at Station M10 are used to initialize the wind 

field for model runs. Table 2 lists the 39 model runs chronologically during these 10 days. The 

surface heating (Q) in table 2 is calculated from equation 24a. The wind speed and direction in 

table 2 are the initialized values which are the 5 min averages of the observed wind at Station 

M10, located in the middle of the diffusion study area. 
18

Table 2. The 39 cases for the HRW simulations. Q is defined by equation 24a in the text. 

CaseNo. Day 

(1992) 

Local 

Time 

Q

(K) 

Wind Speed 

(ms–1) 

WindDirection 

(degree) 

1 14 Sept 14:10' 9.21 5.4 261 

2 15:55' 9.78 7.2 280 

3 17:20' 2.58 4.0 270 

4 15 Sept 12:10' –5.10 8.0 240 

5 13:55' –2.50 7.5 240 

6 15:55' –2.70 7.4 245 

7 17:40' 2.98 5.5 240 

8 16 Sept 10:30' –2.45 3.0 120 

9 13:55' –3.45 0.3 90 

10 16:10' –5.45 3.0 160 

11 17:50' –5.52 2.0 160 

12 17 Sept 10:30' –6.54 5.0 130 

13 14:15' –5.65 1.0 180 

14 17:25' –4.02 5.5 180 

15 18 Sept. 10:30' –0.54 4.5 280 

16 13:10' –0.88 4.0 255 

17 14:55' –2.34 3.0 250 

18 18:15' 0.00 2.5 170 

19 19 Sept 10:30' 3.88 4.0 260 

20 12:10' 6.12 5.0 270 

21 14:00' 6.59 5.0 270 

22 16:10' 5.94 4.5 275 

23 18:25' 3.56 1.5 250 

24 20 Sept 10:30' 3.73 4.5 220 

25 12:40' 10.26 2.5 225 

26 14:45' 5.39 2.5 225 

27 15:55' 3.86 2.0 200 

28 21 Sept 10.30' –0.20 7.5 70 

29 14:25' 1.29 4.5 110 

30 15:25' –1.35 3.5 140 

31 18:00' –1.83 3.5 160 

32 22 Sept 10:30' 1.11 0.5 320 

33 13:30' 7.32 1.2 295 

34 15:10' 2.20 2.5 360 

35 23 Sept 10:30' 5.44 6.0 280 

36 12:00' 8.78 7.5 275 

37 14:15' 8.76 8.5 270 

38 16:10' 9.71 5.0 280 

39 17:40' 1.48 2.5 250 
19

As mentioned before, the domain of the model simulation is 5 by 5 km with a spatial resolution 

of 100 m. Thus, there are 2601 (51 times 51) grid points covered by the computation domain. 

Only the simulated wind speeds and directions at the 11 gridpoints corresponding to the 11 

measurement locations have been chosen for comparison. To evaluate the model performance, 

the simple linear regression analysis is used as follows 

y = ax + b (42) 

where x is the observed value of either wind speed or wind direction, and y is the corresponding 

simulated value from the model runs. The terms a and b represent the slope and intercept of the 

regression line, equation 42, respectively. For each run (case) listed in table 2, there are 

generally 11 pairs (datapoints) of (x, y). However, there are only 10 datapoints for case 9, and 9 

datapoints for case 23. In these two cases, the observed wind speeds from M12 (case 9) and 

from both M12 and M13 (case 23) were zero. Consequently, the wind direction is undetermined 

for these two cases. Therefore, the total number of (x, y) datapoints from the 39 cases is 426 for 

the simultaneous wind direction and speed evaluation. For a single station, except M12 and 

M13, on the other hand, there are 39 datapoints for the linear regression analysis. A correlation 

coefficient (r) and a standard deviation (S) for a set of (x, y) datapoints can also be calculated. 

The linear regression should yield a = 1, b = 0, r = 1, and S = 0 for an ideal or perfect model. In 

reality, however, no model is perfect. The values of a, b, r, and S can help evaluate the model 

performance under various circumstances. 

3.3 The Results of Linear Regression Analysis 

(1) Figures 4a and 4b present the simulated winds from the 39 cases versus the observed winds 

in scatter diagrams with the regression line equation 42 superimposed. Figure 4a is for the wind 

direction comparison, while figure 4b is for the wind speed comparison. Both figures show the 

plot for all 11 stations together as well as for individual stations. Table 3 lists the values of a, b, 

r, and S corresponding to figures 4a and 4b. From the total 426 datapoints, the values of a, b, r, 

and S are 1.004, 2.000, 0.970, and 16.3°, respectively for the wind direction, and 0.868, 0.902, 

0.857, and 1.02 ms–1, respectively for the wind speed, as shown in table 3. This result has been 

partially reported by Cionco and Byers (1995). Therefore, the HRW model simulations as a 

whole appear to provide quite high positive correlations, especially for the wind direction. For 

individual stations, r is always greater than 0.96 for wind direction and greater than 0.80 for 

wind speed, as indicated in table 3. The model performs fairly well for each of the 11 locations, 

although it appears to perform better for certain locations (e.g., M10, M11) than for other 

locations (e.g., M15, M14), as seen from figure 4 and table 3. The model does not demonstrate 

significant failure at any particular location. It is also evident that the model simulations appear 

better for wind direction than for wind speed. As for an individual case, the model sometimes 

simulates a case satisfactorily, but not always. A statistical result for each individual case is 

meaningless due to small number of data points, and hence is not presented. 
20

Figure 4a. Scatter diagrams for simulated wind direction (y, degree) versus observed values (x, degree). 

The straight lines result from the linear regression analysis, equation 42. 
21

Figure 4b. Scatter diagrams for simulated wind speed (y, ms–1) versus observed values (x, ms–1). 

The straight lines result from the linear regression analysis, equation 42. 
22

Table 3. The value of n (data pair number), a (slope), b (intercept), r (correlation coefficient), and S (standard 

deviation) from the linear regression equation 42, for wind direction and wind speed, respectively. 

Wind direction (degree) Wind speed (ms–1 Station n ) 

a b r S a b r S 

M2 39 1.0116 –5.3141 0.9716 15.9 0.8169 0.9626 0.8739 1.02 

M4 39 1.0288 –7.0611 0.9762 15.1 0.7601 1.0301 0.8358 1.00 

M5 39 1.0241 –8.4746 0.9759 14.7 0.7964 0.8659 0.9397 0.57 

M6 39 1.0795 –11.2393 0.9676 17.9 0.8728 0.8877 0.8582 1.07 

M7 39 1.0278 5.9955 0.9697 17.0 0.8448 1.2949 0.8942 0.92 

M10 39 1.0386 –10.035 0.9926 8.3 0.8951 0.5111 0.9586 0.52 

M11 39 0.9833 7.5243 0.9862 11.0 0.3993 0.9655 0.9255 0.73 

M12 37 1.0344 –5.6589 0.9778 14.8 0.8725 0.6401 0.8629 0.90 

M13 38 0.9391 20.405 0.9742 14.9 1.0292 1.0245 0.8502 1.07 

M14 39 0.9631 15.7441 0.9775 14.5 0.8834 0.8918 0.8038 1.40 

M15 39 0.9255 7.0691 0.9613 18.7 1.2376 –0.5478 0.8708 1.11 

All 426 1.0043 2.0001 0.9701 16.3 0.8681 0.9024 0.8567 1.02 

(2) To further investigate the model performance, comparisons of wind field between 

simulations and observations for different wind directions and different wind speeds have been 

conducted. The 39 cases listed in table 2 have been grouped into three categories of wind 

direction and three categories of wind speed. The criteria for wind direction categorization take 

into account the general ridge direction in the MADONA area, which is approximately 230°–

50°, as mentioned earlier (see figure 3). The three wind direction (WD) categories are defined as 

1) Parallel to the ridge (14 cases), 

 200°  WD  260° or 20°  WD  80° (43a) 

2) Perpendicular to the ridge (10 cases) 

 110°  WD  170° or 290°  WD  350° (43b) 

3) Slant to the ridge (15 cases) 

 350°  WD  20° or 80°  WD  110° 

 or 170°  WD  200° or 260°  WD  290° (43c) 

The three wind speed (WS) categories are defined more or less subjectively. They are 

1) Light wind (15 cases), 

 WS ≤ 3.0 ms–1 (43d) 

2) Moderate wind (16 cases), 

 3.0 ms–1  WS  6.0 ms–1 (43e) 
23

3) Strong wind (8 cases), 

 WS ≥ 6.0 ms–1 (43f) 

The results of these comparisons are presented in table 4. They are also shown in figure 5 for the 

three wind direction categories, and in figure 6 for the three wind speed categories, respectively. 

Table 4. Same as table 3 except for three wind direction categories as defined by equation 43 (a, b, c) and three 

wind speed categories as defined by equation 43 (d, e, f). 

Wind Direction 

(degrees) 

Wind Speed 

(ms–1 Category n ) 

a b r S a b r S 

Wind direction: 

Parallel 

Perpendicular 

Slant 

152 

110 

164 

0.9684 

1.0028 

0.9530 

11.7935 

–1.5884 

16.7129 

0.9335 

0.9737 

0.9484 

17.8 

15.4 

14.9 

0.7946 

0.7447 

0.9111 

0.9584 

1.1904 

1.0255 

0.8878 

0.7117 

0.8144 

0.90 

0.83 

1.12 

Wind speed: 

Light 

Moderate 

Strong 

162 

176 

88 

0.9787 

1.0271 

1.0013 

6.5416 

–1.5462 

1.1356 

0.9452 

0.9829 

0.9888 

21.8 

12.4 

10.0 

0.7949 

0.8645 

0.1479 

0.9477 

1.0031 

6.3077 

0.6617 

0.7293 

0.1507 

0.79 

1.04 

1.04 
24

Figure 5. Same as figure 4a except for three wind direction categories as defined by equation 43a, b, c. 
25

Figure 6. Same as figure 4b except for three-wind speed categories as defined, by equation 43d, e, f. 
26

It appears from table 4 and figure 5 that the model performed better for the perpendicular wind 

direction category with the highest value of r (0.974) for wind direction and the lowest value of S

(0.83 ms–1) for wind speed. It is conceivable that the airflow will experience more acceleration or 

deceleration when its direction is perpendicular to a ridge. Consequently, the magnitude related 

to the A



 term in equation 9 will be more significant as compared to the parallel or slant wind 

directions. For the three categories based on wind speed, figure 6 and table 4 reveal an 

unsatisfactory situation for the strong wind category. The scatter diagram at the right lower 

corner of figure 6 shows that the datapoints are scattered around the regression line with low 

value of r (0.1507) and high value of S (1.04 ms–1) although the simulations for the wind 

direction appear much better than for the wind speed. This result receives more discussion in the 

following subsection. 

(3) Finally, the 39 cases have also been grouped in terms of atmospheric stability for model 

performance evaluation. Based on the value of Q, equation 24a, the following three stability 

categories are defined as 

1) Unstable (19 cases),

Q > 2.0° (44a) 

2) Near-neutral (9 cases),

Q  2.0° (44b) 

3) Stable (11 cases),

Q < –2.0° (44c) 

The results for the three stability categories are shownin figure 7 and table 5. Among the three 

atmospheric stability categories, the model appears to perform the best for the near-neutral 

stratification as compared to unstable or stable stratification. For example, both values of r for 

wind direction and speed are the highest for the near-neutral category. This result will be 

discussed in the next section. 
27

Figure 7. Same as figure 4 except for three stability categories as defined by equation 44a, b, c.
28

Table 5. Same as table 3 except for three atmospheric stability categories defined by equation 44a, b, c. 

Category n Wind direction (degree) Wind speed (ms–1) 

 a b r S a b r S 

Unstable 207 0.8949 33.8784 0.8881 15.7 0.9291 0.9063 0.8561 1.00 

Near-neutral 99 1.0112 –0.4089 0.9759 16.4 0.7165 0.8671 0.8711 0.67 

Stable 120 0.9379 11.1184 0.9444 16.0 0.7445 1.4706 0.8512 1.01 

3.4 Discussion 

(1) As presented in the last subsection, the HRW model simulates the wind field better for the 

near-neutral condition than for either stable or unstable conditions. This is because under neutral 

conditions (b = 0) the buoyancy acceleration plays no role and the wind field is influenced by 

topography only. This result implies that the buoyancy (b) term is important and that a realistic 

representation of this term for the model is a challenge. As defined by equations 24 and 24a, the 

calculation of the buoyancy factor (b) or the surface heating (Q) involves the background 

(reference) potential temperature (o), the surface air temperature (Tsfc), and the related surface 

air potential temperature (sfc). The calculation of these three input parameters is not trivial. 

Currently o is calculated through a linear extrapolation from the potential temperature at 700 mb 

and 850 mb obtained from a nearby radiosonde observation. In addition, the current algorithm 

uses Q as input parameter for calculation of Tsfc or sfc , which can be quite variable over 

complex terrain. The HRW model employs a low level of sophistication and does not invoke a 

surface energy balance scheme to deal with surface temperature. It has also been noted that the 

model is essentially a single surface layer model which decouples this layer from the rest of the 

boundary layer. This simplification imposes inherent difficulty for the representation of o.

Lanicci (1985) had suggested a modified technique for calculation of o, Tsfc, and sfc. His 

technique, however, was not verified by observations. Preliminary tests of model sensitivity to 

atmospheric stabilities (changes in Q or b) have been conducted by one of the authors (Huynh). 

Those tests show that the HRW simulation is more sensitive to the change in surface heating or 

surface temperature under stable conditions than it is under unstable conditions. Results of those 

sensitivity studies will be reported elsewhere. 

(2) The results presented in the previous subsection indicate that the HRW model simulates the 

wind direction field with less error than the wind speed field for the MADONA. A brief 

discussion of the relative error of the wind direction (ED) and of the wind speed (ES) appears 

useful. These relative errors and its ratio (RDS) are defined as 

,, , D

D S DS

S

dD dS E E ER

DS E   (45) 
29

where dD and dS denote the errors of wind direction (D) and wind speed (S), respectively, and 

S2

 = u

2

 + v

2

 , D = const. –  ,  = arctan(v/u) , (46) 

where u and v are defined as eastward and northward component of the wind vector, respectively 

according to meteorological convention. Similarly, the relative errors in u and v, and their ratio 

(Ruv) are defined as 

,, . u

u v uv

v

du dv E EER

uvE   (47) 

A relationship between RDS and Ruv can be derived. Equations 48 and 49 can be derived from the 

definitions of equations 45, 46, and 47 

2 2 22

2 2 22 . S uv

u du v dv u v E EE

Su Sv S S

         (48) 

  2

1

arctan . D uv

v uv E d EE

D u DS

               (49) 

Equation 50 may be derived from the two equations above 

   

  2 2 2

1 1 1 tan . tan

u v uv

DS

u v uv

E E uv R

R

DuE vE D R





     

 (50) 

The wind direction (D) in equation 50 is defined in the range: 0 < D < 2. For the MADONA 

simulations, D > /3 = 1.0472, see figure 4. Therefore, the following relation is true 

 

2

1 tan . tan

uv

DS

uv

R

R

R





 



 (51) 

Note that Ruv ≥ 0 if du and dv result only from wind speed errors. A positive error in wind speed 

(dS > 0) causes positive Eu and Ev while a negative error (dS < 0) results in negative Eu and Ev. 

Consequently, the absolute value of RDS can be estimated from equation 51. 

-1

uv

1 for tan 1

1-R for tan 1

uv

DS

R

R





 





 (52) 

The above expression implies that RDS is <1.0 when 0.5 < Ruv < 2.0. Only if Ruv is either large 

or small, out of the above range, can the relative errors in wind direction estimation be greater 

than the relative errors in wind speed estimation. This situation is likely to occur when the wind 

direction is east, west (v = 0), south, or north (u = 0). For the MADONA simulations, such wind 

directions occur infrequently. Consequently, the statistical results for wind direction simulations 

appear better than for wind speed simulations. The above argument, however, is not intended to 
30

provide a physical explanation for our results. Cauchey et al. (1979), for example, have 

formulated a wind direction budget equation. A similar wind speed budget equation is not 

difficult to derive. It is, however, beyond the scope of the present study to analyze these physical 

budgets. 

(3) Finally, the performance of the model simulations for the strong wind speed category seems 

undesirable. There are eight cases in this category, equation 43c. They are cases 2, 4, 5, 6, 28, 

35, 36, and 37 (table 2). As shown in figure 6 (the lowest right panel), the datapoints are 

scattered around the regression line. The correlation coefficient (r) for the wind speed is 0.1507 

(table 4). The difference between the simulated and observed wind can be as large as 27 degrees 

in wind direction and 3.5 ms–1 in wind speed, the extremes at an individual station. 

A careful examination of each of these eight cases indicates that the wind fields at the time of the

simulations all exhibit significant unsteadiness. There were considerable variations 

(accelerations) of wind vector, especially wind speed, around the time of the simulations, i.e., 

considerable   V t /

 . If there is a gusty wind under a strong wind speed situation, a 1 ms–1 of 

wind speed change in a period of 5 min, then this fluctuation alone is equivalent to an 

acceleration on the order of 10–3 ms–2, which can be greater than the magnitude of acceleration 

induced by topography or buoyancy. Recall that the HRW model is not a prognostic model, just 

a diagnostic one, which is applied to a steady wind field. It comes with no surprise then that the 

model will perform less satisfactorily for gusty (unsteady) wind fields. 

4. Conclusions 

We have presented a comprehensive description of the high resolution wind model which has 

been an undertaking at the U.S. Army Research Laboratory (formerly U.S. Army Atmospheric 

Sciences Laboratory) for three decades. From the early effort of model development, several 

requirements have been kept in mind: (1) the model can be applied over very complex terrain 

with limited meteorological information (inputs); (2) the model is adaptable to small horizontal 

grid spacing on the order of 100 m for a horizontal domain smaller than 10 by 10 km; (3) the 

model is to be run as a combined system with existing atmospheric dispersion and transport 

models; (4) the model can be used on a small computer. Our HRW model uses Gauss’s Principle 

of least constraint for a variational adjustment of an initial estimated wind field in a single 

surface layer to conform with terrain structures, mass conservation, and buoyancy forces. This 

basic approach and its features appear to meet the above requirements. 

The measurements from the field study of Meteorology and Diffusion Over Non-Uniform Areas 

(MADONA) have provided very valuable data to test our HRW model. As shown in the last 

section, the results from the present study demonstrate that the HRW model generally performs 

well for a total of 39 cases when compared to MADONA data. The results of our evaluation, on 
31

the other hand, also indicate certain limitations in the HRW model applications. For example, 

the model may not simulate a gusty (unsteady) wind field as successfully as a steady wind field. 

Examination of these limitations and further testing of the model will allow us to gain insights 

consistent with the model physics and to develop a next generation model. 

Finally, the results of a companion evaluation of the HRW model (Williamson et al. 2005) 

indicate that the correlation coefficient by itself is only a weak indicator of the model validity. In 

fact, that study showed that the HRW model performance for the MADONA field data is 

actually slightly worse than the homogeneous wind field used to initiate it, and this holds true for 

all statistical measures that were considered. 
32

5. References 

Arya, S. P. Air Pollution Meteorology and Dispersion, Oxford University Press, 1999, 98–100. 

Ball, J. A.; Johnson, S. A. Physically-based High-resolution Surface Wind and Temperature 

Analysis for EPAMS; ASL-CR-78-0043-1; 1978, 270. 

Bradshaw, P. The Analogy Between Streamline Curvature and Buoyancy in Turbulent Shear 

Flow. J. Fluid Mech. 1969, 36, Part 1,177–191. 

Businger, J. A.; Wyngaard, J. C.; Izumi, Y.; Bradley, E. F. Flux Profile Relationships in the 

Atmospheric Boundary Layer. J. Atmos. Sci. 1971, 28, 181–189. 

Caughey, S. J.; Wyngaard, J. C.; Kaimal, J. C. Turbulence in the Evolving Stable Boundary 

Layer. J. Atmos. Sci. 1979, 36, 1041–1052. 

Cionco, R. M. On the Coupling of Canopy Flow to Ambient Flow for a Variety of Vegetation 

Types and Densities. Boundary-Layer Meteorology 1983, 26, 325–335. 

Cionco, R. M. Modeling Windfields and Surface Layer Profiles Over Complex Terrain and 

Within Vegetative Canopies, in the Forest-Atmosphere Interactions; Eds. B. A. Hutchinson 

and B. B. Hicks, D. Reidel Publishing Company, 1985. 

Cionco, R. M. Windfield Simulations from a Hierarchy of Nested Meso- and 

Micrometeorological Models. Proc. Sympo. Mesoscale Analysis & Forecasting, Vancouver, 

Canada, 17–19 August 1987. 

Cionco, R. M.; Byers, J. H. Initialising the HRW Model with Meteorological Inputs from 

Different Stes to Optimize Model Output, Preprint Volume of the 11th Symposium on 

Boundary Layers and Turbulence, 27–31 March 1995, Charlotte, NC by the American 

Meteorological Society, Boston, MA, 1995. 

Cionco, R. M.; Byers, J. H. Influence of Airfield Structures and Complex Terrain upon 

Simulated Wind Fields, Preprint Volume of the 12th Symposium on Boundary Layers and 

Tubulence, 28 July–1 August 1997, Vancouver, BC Canada by American Meteorological 

Society, Boston MA, 1997.

Cionco, R. M. Simulating Wind Fields for Different Resolutions of Terrain and Morphological 

Land Features, Second Urban Environment Symposium, 13th Conference on 

Biometeorology and Aerobiology, 2–6 November 1998, Albuquerque, New Mexico, 

American Meteorological Society, 1998. 
33

Cionco, R. M.; Weber, H;. aufm Kampe, W. High Resolution Meteorological Transport and 

Diffusion Simulations in a Complex river Valley, Preprint volume of the 23rd Conf. on 

Agricultural and Forest Meteorology, 13th conf. on Biometeorology and Aerobiology, and 

2nd Urban Environment Symp., 2–6 November 1998, Albuquerque, NM, by AMS, Boston, 

MA, 1998. 

Cionco, R. M.; Ellefsen, R. High Resolution Urban Morphology Data for Urban Wind Flow 

Modeling. Atmospheric Environment 1998, 32 (1), 7–17. 

Cionco, R. M.; et al. An Overview of MADONA: A Multinational Field Study of HighResolution Meteorology and Diffusion over Complex Terrain. Bulletin of the American 

Meteorological Society 1999, 80 (1), 5–19. 

Garratt, J. R. The Atmospheric Boundary Layer, Cambridge University Press, 1992, 316. 

Lanczos, C. The Variational Principles of Mechanics (Third Edition), University of Toronto 

Press, 1966, 106–110. 

Lanicci, J. M., Sensitivity Tests of a Surface-Layer Windflow Model to effects of stability and 

Vegetation; AFGL-TR-85-0265; Air Force Geophysics Laboratory: Hanscom AFB, MA 

01731, 1985. 

Orlanski, I. A rational Subdivision of Scales of Atmospheric Processes. Bulletin of the Amer. 

Meteor. Soc. 1975, 56, 527–530. 

Panofsky, H. A.; Dutton, J. A. Atmospheric Turbulence; John Wiley & Sons, Inc. 1984; 397. 

Stull, R. B. An Introduction to Boundary Layer Meteorology; Kluwer Academic Publishers, 

1988; 666. 

Tai, Chen-To. Generalized Vector and Dyadic Analysis. IEEE Press 1992, 134. 

Thykier-Nielsen S.; Mikkelsen, T.; Santabarbara, J.; Ride, D. J.; Higgs, T.; Weber, H. 

MADONA: Real-Time Diffusion Model Simulations, Preprint volume of the 11th Symp. on 

Boundary Layers and Turbulence, 27–31 March 1995, Charlotte, NC by AMS, Boston, MA, 

1995. 

Wallace, J. M.; Hobbs, P. V. Atmospheric Science, An Introductory Survey; Academic Press, 

1977; 467. 

Weber, H.; aufm Kampe, W.; Cionco, R. M. Visualization of the MADONA Data Base and Use 

of Selected Sequences in a Wind Flow and Diffusion Simulation System, Preprint volume of 

the 11th Symp. on Boundary Layers and Turbulence, 27–31 March 1995, Charlotte, NC, by 

AMS, Boston, MA, 1995. 
34

Williamson, C. C.; Garvey, D. M.; Chang, S. S. Revisiting the Analysis of the High Resolution 

Wind Model Using MADONA Field Data; ARL-TR-3497; U.S. Army Research Laboratory: 

Adelphi, MD, 2005. 
34

No. of

Copies Organization

1 ADMNSTR 

ELECT DEFNS TECHL INFO CTR 

 ATTN DTIC OCP 

 8725 JOHN J KINGMAN RD STE 0944 

 FT BELVOIR VA 22060-6218 

1 CHAIRMAN JOINT CHIEFS OF 

 STAFF 

 ATTN J5 R&D DIV 

 WASHINGTON DC 20301 

1 DARPA 

 ATTN IXO S WELBY 

 3701 N FAIRFAX DR 

 ARLINGTON VA 22203-1714 

2 DIR OF DEFNS RSRCH & ENGRG 

 ATTN DD TWP 

 ATTN ENGRG 

 WASHINGTON DC 20301 

1 CD OFC OF THE SECY OF DEFNS 

 ATTN ODDRE (R&AT) 

 THE PENTAGON 

 WASHINGTON DC 20301-3080 

1 COMMANDING OFFICER 

 ATTN NMCB23 

 6205 STUART RD STE 101 

 FT BELVOIR VA 22060-5275 

1 US ARMY RSRCH DEV AND ENGRG 

 CMND 

 ARMAMENT RSRCH DEV AND 

 ENGRG CTR 

 ARMAMENT ENGRG AND 

 TECHNLGY CTR 

 ATTN AMSRD AAR AEF T 

 J MATTS 

 BLDG 305 

 ABERDEEN PROVING GROUND MD 

 21005-5001 

1 DIR OF CHEM & NUC OPS DA 

 DCSOPS 

 ATTN TECHL LIB 

 WASHINGTON DC 20301 

1 NATL GROUND INTLLGNC CTR 

 ATTN RSRCH & DATA BRANCH 

 CHARLOTTESVILLE VA 22901-5396 

No. of 

Copies Organization

1 PM TIMS, PROFILER (MMS-P) 

 AN/TMQ -52 

 ATTN B GRIFFIES 

 BUILDING 563 

 FT MONMOUTH NJ 07703 

1 US ARMY ENGRG DIV 

 ATTN HNDED FD 

 PO BOX 1500 

 HUNTSVILLE AL 35807 

1 US ARMY INFO SYS ENGRG CMND 

 ATTN AMSEL IE TD F JENIA 

 FT HUACHUCA AZ 85613-5300 

1 US ARMY MIS & SPC INTLLGNC CTR 

 AIAMS YDL 

 REDSTONE ARSENAL AL 35898-5500 

1 US ARMY NUC & CHEML AGCY 

 7150 HELLER LOOP STE 101 

 SPRINGFIELD VA 22150-3198 

1 COMMANDER 

 US ARMY RDECOM 

 ATTN AMSRD AMR W C MCCORKLE 

 5400 FOWLER RD 

 REDSTONE ARSENAL AL 35898-5000 

1 US ARMY STRTGC DEFNS CMND 

 ATTN CSSD H MPL TECHL LIB 

 PO BOX 1500 

 HUNTSVILLE AL 35807 

1 CHIEF OF NAV OPS DEPT OF THE 

 NAVY 

 ATTN OP 03EG 

 WASHINGTON DC 20350 

1 US AIR FORCE TECH APPL CTR 

 ATTN HQ AFTAC/TCC 

 1030 SOUTH HIGHWAY A1A 

 PATRICK AFB FL 32925-3002 

1 CENTRAL INTLLGNC AGCY DIR DB 

 STANDARD 

 ATTN OSS/KPG/DHRT 

 1E15 OHB 

 WASHINGTON DC 20505 
35

No. of

Copies Organization

1 US GOVERNMENT PRINT OFF 

 DEPOSITORY RECEIVING SECTION 

 ATTN MAIL STOP IDAD J TATE 

 732 NORTH CAPITOL ST NW 

 WASHINGTON DC 20402 

1 NATL CTR FOR ATMOS RSRCH 

 ATTN NCAR LIBRARY SERIALS 

 PO BOX 3000 

 BOULDER CO 80307-3000 

1 US ARMY RSRCH LAB 

 ATTN AMSRD ARL CI OK TP TECH 

 LIB T LANDFRIED 

 BLDG 4600 

 ABERDEEN PROV 21005-5066 

1 DIRECTOR 

 US ARMY RSRCH LAB 

 ATTN AMSRD ARL RO EV 

 W D BACH 

 PO BOX 12211 

 RESEARCH TRIANGLE PARK NC 

 27709 

7 US ARMY RSRCH LAB 

 ATTN AMSRD ARL CI ED 

 D HOOCK (5 COPIES) 

 ATTN AMSRD ARL CI D KNAPP 

 (2 COPIES) 

 BLDG 1622 

 WHITE SANDS MISSILE RANGE NM 

 88002-5501 

No. of

Copies Organization

39 US ARMY RSRCH LAB 

 ATTN AMSRD ARL CI E P CLARK 

 ATTN AMSRD ARL CI ED 

 C WILLIAMSON (5 COPIES) 

 ATTN AMSRD ARL CI ED 

 D GARVEY 

 (5 COPIES) 

 ATTN AMSRD ARL CI ED 

 G D HUYNH (5 COPIES) 

 ATTN AMSRD ARL CI ED 

 S CHANG (20 COPIES) 

 ATTN AMSRD ARL CI OK PE 

 TECHL PUB 

 ATTN AMSRD ARL CI OK TL 

 TECH LIB 

 ATTN IMNE ALC HRR MAIL & 

 RECORDS MGMT 

 ADELPHI MD 20783-1197 

TOTAL: 70 (68 HCs, 1 CD, 1 ELECT) 
36

INTENTIONALLY LEFT BLANK