Technotes 576
Technotes 576
Advanced Research
WRF Model Version 4
William C. Skamarock
Joseph B. Klemp
Jimy Dudhia
David O. Gill
Zhiquan Liu
Judith Berner
Wei Wang
Jordan G. Powers
Michael G. Duda
Dale M. Barker
Xiang-Yu Huang
NCAR/TN-556+STR
National Center for
Atmospheric Research
P. O. Box 3000
Boulder, Colorado
80307-3000
www.ucar.edu
NCAR TECHNICAL NOTES
http://library.ucar.edu/research/publish-technote
IA – Instructional Aids
Instruction manuals, bibliographies, film supplements,
and other research or instructional aids.
PROC – Proceedings
Documentation or symposia, colloquia, conferences,
workshops, and lectures. (Distribution maybe limited to
attendees).
A Description of the
Advanced Research WRF Model Version 4
William C. Skamarock
Joseph B. Klemp
Jimy Dudhia
David O. Gill
Zhiquan Liu
Judith Berner
Wei Wang
Jordan G. Powers
Michael G. Duda
Dale M. Barker
Xiang-Yu Huang
Preface xi
Acknowledgments xiii
1 Introduction 1
1.1 Advanced Research WRF (ARW) . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Major Features of the ARW System, Version 4 . . . . . . . . . . . . . . . . . . . 3
2 Governing Equations 7
2.1 Vertical Coordinate and Flux-Form Variables . . . . . . . . . . . . . . . . . . . . 7
2.2 Flux-Form Euler Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3 Map Projections, Coriolis and Curvature Terms . . . . . . . . . . . . . . . . . . 10
2.4 Perturbation Form of the Governing Equations . . . . . . . . . . . . . . . . . . . 12
3 Model Discretization 15
3.1 Temporal Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.1.1 Runge-Kutta Time Integration Scheme . . . . . . . . . . . . . . . . . . . 15
3.1.2 Acoustic Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.1.3 Full Time-Split Integration Sequence . . . . . . . . . . . . . . . . . . . . 18
3.1.4 Diabatic Forcing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.1.5 Hydrostatic Option . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.2 Spatial Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.2.1 Acoustic Step Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.2.2 Coriolis and Curvature Terms . . . . . . . . . . . . . . . . . . . . . . . . 22
3.2.3 Advection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.2.4 Pole Conditions for the Global Latitude-Longitude Grid . . . . . . . . . 27
3.3 Stability Constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.3.1 RK3 Time Step Constraint . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.3.2 Acoustic Time Step Constraint . . . . . . . . . . . . . . . . . . . . . . . 28
3.3.3 Adaptive Time Step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.3.4 Map Projection Considerations . . . . . . . . . . . . . . . . . . . . . . . 30
i
4.2.2 Horizontal and Vertical Diffusion in Physical Space . . . . . . . . . . . . 33
4.2.3 Computation of the Eddy Viscosities . . . . . . . . . . . . . . . . . . . . 35
4.2.4 TKE equation for the 1.5 Order Turbulence Closure . . . . . . . . . . . . 36
4.2.5 Sixth-Order Spatial Filter on Coordinate Surfaces . . . . . . . . . . . . . 37
4.3 Filters for the Time-split RK3 scheme . . . . . . . . . . . . . . . . . . . . . . . . 38
4.3.1 Three-Dimensional Divergence Damping . . . . . . . . . . . . . . . . . . 38
4.3.2 External Mode Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4.3.3 Semi-Implicit Acoustic Step Off-centering . . . . . . . . . . . . . . . . . . 39
4.4 Formulations for Gravity Wave Absorbing Layers . . . . . . . . . . . . . . . . . 39
4.4.1 Absorbing Layer Using Spatial Filtering . . . . . . . . . . . . . . . . . . 40
4.4.2 Implicit Rayleigh Damping for the Vertical Velocity . . . . . . . . . . . . 40
4.4.3 Traditional Rayleigh Damping Layer . . . . . . . . . . . . . . . . . . . . 41
4.5 Other Damping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.5.1 Vertical-Velocity Damping . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5 Initial Conditions 43
5.1 Initialization for Idealized Conditions . . . . . . . . . . . . . . . . . . . . . . . . 43
5.2 Initialization for Real-Data Conditions . . . . . . . . . . . . . . . . . . . . . . . 45
5.2.1 Use of the WRF Preprocessing System by ARW . . . . . . . . . . . . . . 45
5.2.2 Reference State . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
5.2.3 Vertical Interpolation and Extrapolation . . . . . . . . . . . . . . . . . . 47
5.2.4 Perturbation State . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.2.5 Generating Lateral Boundary Data . . . . . . . . . . . . . . . . . . . . . 48
5.2.6 Masking of Surface Fields . . . . . . . . . . . . . . . . . . . . . . . . . . 49
5.3 Digital Filtering Initialization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.3.1 Filter Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.3.2 DFI Schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
5.3.3 Backward Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
7 Nesting 57
7.1 Horizontal Nesting Options . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
7.2 Staggering and Feedback . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
7.3 Nested Lateral Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . 64
7.4 Steps to Generate a Nest Grid . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
8 Physics 69
8.1 Microphysics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
8.1.1 Kessler scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
8.1.2 Purdue Lin scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
ii
8.1.3 WRF Single-Moment 3-class (WSM3) scheme . . . . . . . . . . . . . . . 70
8.1.4 WSM5 scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
8.1.5 WSM6 scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
8.1.6 WDM5 and WDM6 schemes . . . . . . . . . . . . . . . . . . . . . . . . . 71
8.1.7 WSM7 and WDM7 schemes . . . . . . . . . . . . . . . . . . . . . . . . . 72
8.1.8 Eta Grid-scale Cloud and Precipitation (2001) scheme . . . . . . . . . . . 72
8.1.9 Thompson et al. and aerosol-aware Thompson-Eidhammer schemes . . . 72
8.1.10 Goddard Cumulus Ensemble Model scheme . . . . . . . . . . . . . . . . 72
8.1.11 Goddard 4ICE scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
8.1.12 Morrison et al. 2-Moment scheme . . . . . . . . . . . . . . . . . . . . . . 73
8.1.13 Milbrandt-Yau Double-Moment scheme . . . . . . . . . . . . . . . . . . . 74
8.1.14 CAM Morrison-Gettelman scheme . . . . . . . . . . . . . . . . . . . . . . 74
8.1.15 Stony-Brook University Lin-Colle scheme . . . . . . . . . . . . . . . . . . 74
8.1.16 NSSL microphysics schemes . . . . . . . . . . . . . . . . . . . . . . . . . 74
8.1.17 HUJI spectral bin microphysics schemes . . . . . . . . . . . . . . . . . . 74
8.1.18 Predicted Particle Properties (P3) scheme . . . . . . . . . . . . . . . . . 74
8.1.19 Jensen ISHMAEL microphysics . . . . . . . . . . . . . . . . . . . . . . . 75
8.2 Cumulus Parameterization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
8.2.1 Kain-Fritsch schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
8.2.2 Betts-Miller-Janjic scheme . . . . . . . . . . . . . . . . . . . . . . . . . . 76
8.2.3 Grell Schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
8.2.4 Simplified Arakawa-Schubert Schemes . . . . . . . . . . . . . . . . . . . . 78
8.2.5 Zhang-McFarlane Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . 79
8.2.6 Tiedtke Schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
8.3 Shallow Cumulus Parameterization . . . . . . . . . . . . . . . . . . . . . . . . . 80
8.3.1 University of Washington Scheme . . . . . . . . . . . . . . . . . . . . . . 80
8.3.2 GRIMs Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
8.3.3 NSAS Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
8.3.4 Deng Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
8.4 Surface Layer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
8.4.1 Revised MM5 similarity theory . . . . . . . . . . . . . . . . . . . . . . . 81
8.4.2 Similarity theory (MYJ/Eta) . . . . . . . . . . . . . . . . . . . . . . . . 81
8.4.3 QNSE similarity theory . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
8.4.4 MYNN surface layer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
8.4.5 Similarity theory (PX) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
8.4.6 TEMF surface layer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
8.4.7 Similarity theory (MM5) – old version . . . . . . . . . . . . . . . . . . . 82
8.5 Land-Surface Model and Other Surface Options . . . . . . . . . . . . . . . . . . 82
8.5.1 5-layer thermal diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
8.5.2 Noah LSM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
8.5.3 NoahMP LSM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
8.5.4 Rapid Update Cycle (RUC) Model LSM . . . . . . . . . . . . . . . . . . 83
8.5.5 Pleim-Xiu LSM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
8.5.6 Community Land Model (CLM4) . . . . . . . . . . . . . . . . . . . . . . 84
8.5.7 Simplified Simple Biosphere Model (SSiB) . . . . . . . . . . . . . . . . . 85
iii
8.5.8 Urban Canopy Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
8.5.9 Building Environment Parameterization (BEP) . . . . . . . . . . . . . . 85
8.5.10 Building Energy Model (BEM) . . . . . . . . . . . . . . . . . . . . . . . 85
8.5.11 Ocean Mixed-Layer Model . . . . . . . . . . . . . . . . . . . . . . . . . . 86
8.5.12 3-D Ocean Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
8.5.13 CLM4.5 Lake Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
8.5.14 Sea-Ice Treatment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
8.5.15 Updating Lower Boundary Conditions . . . . . . . . . . . . . . . . . . . 87
8.6 Planetary Boundary Layer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
8.6.1 Yonsei University (YSU) PBL . . . . . . . . . . . . . . . . . . . . . . . . 87
8.6.2 Mellor-Yamada-Janjic (MYJ) PBL . . . . . . . . . . . . . . . . . . . . . 88
8.6.3 Quasi-Normal Scale Elimination (QNSE) scheme with EDMF . . . . . . 88
8.6.4 Mellor-Yamada-Nakanishi-Niino (MYNN) Levels 2.5 and 3 . . . . . . . . 88
8.6.5 Asymmetrical Convective Model version 2 (ACM2) PBL . . . . . . . . . 89
8.6.6 Bougeault-Lacarrere PBL . . . . . . . . . . . . . . . . . . . . . . . . . . 89
8.6.7 University of Washington (UW) PBL . . . . . . . . . . . . . . . . . . . . 89
8.6.8 Total Energy Mass Flux (TEMF) PBL . . . . . . . . . . . . . . . . . . . 89
8.6.9 Shin-Hong PBL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
8.6.10 Grenier-Bretherton-McCaa (GBM) PBL . . . . . . . . . . . . . . . . . . 90
8.6.11 Medium Range Forecast Model (MRF) PBL . . . . . . . . . . . . . . . . 90
8.6.12 Gravity Wave Drag . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
8.7 Atmospheric Radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
8.7.1 Rapid Radiative Transfer Model (RRTM) Longwave . . . . . . . . . . . . 91
8.7.2 CAM3 Longwave and Shortwave . . . . . . . . . . . . . . . . . . . . . . . 91
8.7.3 RRTMG Longwave and Shortwave . . . . . . . . . . . . . . . . . . . . . 91
8.7.4 RRTMG-K Longwave and Shortwave . . . . . . . . . . . . . . . . . . . . 92
8.7.5 Eta Geophysical Fluid Dynamics Laboratory (GFDL) Longwave and Short-
wave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
8.7.6 MM5 (Dudhia) Shortwave . . . . . . . . . . . . . . . . . . . . . . . . . . 92
8.7.7 Old Goddard Shortwave . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
8.7.8 New Goddard Longwave and Shortwave . . . . . . . . . . . . . . . . . . 93
8.7.9 Fu-Liou-Gu (FLG) Longwave and Shortwave . . . . . . . . . . . . . . . . 93
8.8 Physics Interactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
10 Nudging 101
10.1 Grid Nudging or Analysis Nudging . . . . . . . . . . . . . . . . . . . . . . . . . 101
10.2 Surface Analysis Nudging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
10.3 Flux-Adjusting Surface Analysis Nudging . . . . . . . . . . . . . . . . . . . . . . 103
iv
10.4 Observational or Station Nudging . . . . . . . . . . . . . . . . . . . . . . . . . . 103
10.5 Spectral Nudging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
C Acronyms 123
References 126
v
vi
List of Figures
3.1 Time step integration sequence. Here n represents the number of acoustic time
steps for a given substep of the RK3 integration, and ns is the ratio of the RK3
time step to the acoustic time step for the second and third RK3 substeps. . . . 19
3.2 Horizontal and vertical grids of the ARW . . . . . . . . . . . . . . . . . . . . . . 21
3.3 Latitude-longitude grid structure in the pole region. In the ARW formulation,
re ∆ψ = ∆y/my . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
5.1 Schematic showing data flow and program components in WPS, and how WPS
feeds initial data to ARW. Text in the rectangular boxes indicates program names.
GEOGRID: defines the model domain and creates static files of terrestrial data.
UNGRIB: decodes GRIB-formatted data. METGRID: interpolates meteorologi-
cal data to the model domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
5.2 An illustration showing the three available DFI schemes: digital filter launch,
diabatic digital filter initialization, and twice digital filter initialization. . . . . . 51
6.1 Specified and relaxation zones for a grid with a single-specified row and column,
plus four rows and columns for the relaxation zone. These are typical values used
for a specified lateral boundary condition for a real-data case. In this figure, the
weight of the relaxation term would be identically zero for the fifth row or column
in from the boundary edge. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
7.1 1-way (parent provides boundary data to child) and 2-way (parent provides
boundary data to child, then child feeds information back to the parent) nesting
options in ARW. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
7.2 Various nest configurations for multiple grids. (a) Telescoping nests. (b) Nests at
the same level with respect to a parent grid. (c) Overlapping grids: not allowed
with feedback activated. (d) Inner-most grid has more than one parent grid: not
allowed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
vii
7.3 Arakawa-C grid staggering for a portion of a parent domain, with an imbedded
nest domain using a 3:1 grid size ratio. Solid lines denote coarse grid cell bound-
aries, and dashed lines are the boundaries for each fine grid cell. The horizontal
components of velocity (“U” and “V”) are defined along the normal cell face, and
the thermodynamic variables (“θ”) are defined at the center of the grid cell (each
square). The bold typeface variables along the interface, between the coarse and
the fine grid, define the locations where the specified lateral boundaries for the
nest are in effect. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
7.4 Similar to Fig. 7.3, but with a 2:1 grid-distance ratio. . . . . . . . . . . . . . . 63
7.5 Nest grid integration sequence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
7.6 Zones of topographic blending for a fine grid. In the fine grid, the first zone is
entirely interpolated from the coarse-grid topography. In the second zone, the
topography is linearly weighted between the coarse grid and the fine grid. . . . . 66
11.1 Sketch showing the relationship between datasets (circles), and algorithms (rect-
angles) of ARW system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
11.2 Sketch of the role of Stage 0 converters in transforming model-specific data (e.g.,
ARW, KMA global model, etc.) to standard perturbation fields and relevant
metadata (e.g., latitude, height, land/sea, etc.). . . . . . . . . . . . . . . . . . . 115
viii
List of Tables
3.1 Maximum stable Courant numbers for one-dimensional linear advection. From
Wicker and Skamarock (2002). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
5.1 Ideal Cases. Listed are the available idealized cases for the Advanced Research
WRF. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
ix
x
Preface
The Advanced Research WRF (ARW) model is a configuration of the Weather Research
and Forecasting (WRF) model. This technical note describes the scientific and algorithmic ap-
proaches in the ARW Version 4, including its dynamical solver, physics options, initialization
capabilities, boundary conditions, grid-nesting techniques, and data assimilation capabilities.
ARW is supported as a community model by the National Center for Atmospheric Research,
facilitating system development and broad use for research, operations, and education. The mod-
eling system supports atmospheric simulations across scales from large-eddy to global. ARWs
applications include real-time NWP, weather events and atmospheric-process studies, data as-
similation development, parameterized-physics development, regional climate simulation, air
quality modeling, atmosphere-ocean coupling, and idealized- atmosphere studies.
This particular version of the Tech Note covers ARW releases up to Version 4.1. This
document will be updated as new releases become available and new features are added to the
model.
xi
xii
Acknowledgments
The Advanced Research WRF (ARW) Model is a community atmospheric modeling system,
and its development and capabilities are the result of the contributions of a host of individuals
and institutions from over the years. The features in ARW Version 4 stem from this gener-
ous and collaborative approach and offer the community a state-of-the-art numerical weather
prediction model reflecting advancements in meteorology, computational science, and data as-
similation. Many people have worked on and provided features and improvements for this latest
version, particularly in the areas of physics and dynamics, and the authors thank them for their
contributions.
The effort to support and maintain ARW is led by the National Center for Atmospheric
Research (NCAR), which is primarily supported by the National Science Foundation (NSF).
The authors sincerely thank NCAR and NSF for their support of the ARW effort. During
the development of Version 4, support in various forms has also been provided by: the National
Oceanic and Atmospheric Administration, the U.S. Department of Defense, the Federal Aviation
Administration, the U.S. Department of Energy, and the Central Weather Bureau of Taiwan.
Lastly, the authors would like to thank the following individuals for their careful review of
this document: George Bryan, Jonathan J. Guerrette, Kevin Manning, Hugh Morrison, Ned
Patton, Brian Reen, Juanzhen Sun, Kelly Werner, and May Wong.
xiii
xiv
Chapter 1
Introduction
The Weather Research and Forecasting (WRF) Model is an atmospheric modeling system de-
signed for both research and numerical weather prediction. WRF is an open-source community
model, and it has been adopted for research at universities and governmental laboratories, for
operational forecasting by governmental and private entities, and for commercial applications by
industry. WRF development began in the latter half of the 1990’s with the goals being to build
a system shared by research and operations and to create a next-generation numerical weather
prediction (NWP) capability. The new modeling system has become a common platform on
which the broad research community can develop capabilities that can transition to operations,
while the extra scrutiny of performance in operations could guide and accelerate development.
The WRF system was developed through a partnership of the National Center for Atmospheric
Research (NCAR), the National Oceanic and Atmospheric Administration (NOAA) (represented
by the National Centers for Environmental Prediction (NCEP) and the NOAA Earth System
Research Laboratory (ESRL)), the United States Air Force, the Naval Research Laboratory, the
University of Oklahoma, and the Federal Aviation Administration. For more information on
the history of the WRF model development see Powers et al. (2017).
1
Static Data
WRF Software Framework
WRF Digital
Preprocessing Filter
System ARW
Dynamical Post Processors
Core Verification
WRFDA Data
Assimilation
Analyses
WRF-Chem
Forecasts
Physics Interface WRF-Hydro
Observations
WRF-Fire
(others)
Physics Packages
from large-eddy to global. ARW’s applications include real-time NWP, weather events and
atmospheric-process studies, data assimilation development, parameterized-physics develop-
ment, regional climate simulation, air quality modeling, atmosphere-ocean coupling, and idealized-
atmosphere studies.
Figure 1.1 depicts the principal components of the ARW system. The WRF Software Frame-
work (WSF) is the infrastructure that contains the dynamics solver, physics packages, utilities
for initialization, WRFDA, and integrated capabilities such as WRF-Chem, WRF-Hydro, WRF-
Fire, etc. The WSF also contains the NMM-E dynamics solver, which is used by NCEP in the
operational HWRF model. The Mesoscale and Microscale Meteorology Laboratory of NCAR
provides support for the ARW, and oversees the WRF repository and releases.
This technical note focuses on the scientific and algorithmic approaches in the ARW Version
4, including its dynamical solver, physics options, initialization capabilities, boundary condi-
tions, and grid-nesting techniques. The WSF provides the software infrastructure, although this
infrastructure is not reviewed in this technical note. Additionally, while WRF-Chem, WRF-
Hydro, WRF-Fire, and other tailored systems use the ARW solver, they are also not covered by
this technical note. For information on actually running the ARW system, the ARW Version 4
Modeling System User’s Guide covers model operation
The following section highlights the major features of the ARW Version 4, first released in
May 2018.
2
1.2 Major Features of the ARW System, Version 4
Model Physics
• Microphysics: Schemes ranging from simplified physics suitable for idealized studies to mixed-
phase, multi-moment, bin, and aerosol-aware approaches to support process studies and
accurate NWP.
3
• Cumulus parameterizations: Deep and shallow convection, adjustment, mass-flux, and scale-
aware schemes available.
• Surface physics: Multi-layer land surface models ranging from a simple thermal model to full
vegetation and soil moisture models, including snow cover and sea ice. Urban parameter-
izations are available.
• Planetary boundary layer physics: Turbulent kinetic energy prediction or non-local K schemes.
• Atmospheric radiation physics: Longwave and shortwave schemes with multiple spectral
bands and a simple shortwave scheme suitable for climate and weather applications. Cloud
effects and surface fluxes are included.
WRFDA System
4
• Portable across a range of available computing platforms.
• Support for multiple physics modules.
• Separation of scientific codes from parallelization and other architecture-specific aspects.
• Input/output Application Program Interface (API) enabling various external packages to be
installed with WRF, thus allowing WRF to easily support various data formats.
• Efficient execution on a range of computing platforms (distributed and shared memory, vector
and scalar types).
WRF-Chem is a full online atmospheric chemistry model with many options and some in-
teractions with the physics via aerosols affecting radiation and microphysics when appropriate
options are chosen. Typically it requires emission source maps as additional inputs. Further
information on WRF-Chem can be found at https://ruc.noaa.gov/wrf/wrf-chem/
WRF-Hydro is a surface hydrological model that can be run online with WRF or offline,
interacting with WRF through the land-surface model. This model can calculate the land water
budget terms including streamflow given additional data for routing. It operates on a higher-
resolution sub-grid relative to the atmospheric model. Further information on WRF-Hydro can
be found at https://ral.ucar.edu/projects/wrf hydro/overview
WRF-Fire operates as a coupled wildland fire model that keeps track of a fire front at
the sub-grid level. It interacts with the model via winds and heat fluxes and requires fuel
information as an additional input. Further information concerning WRF-Fire can be found
at http://www2.mmm.ucar.edu/wrf/users/fire/wrf-fire.html and in the ARW Users Guide at
http://www2.mmm.ucar.edu/wrf/users/docs/user guide v4/contents.html
Additionally, a brief description of these components, along with references to further infor-
mation, can by found in Powers et al. (2017).
5
6
Chapter 2
Governing Equations
The ARW dynamics solver integrates the compressible, nonhydrostatic Euler equations. The
equations are cast in flux form using variables that have conservation properties, following the
philosophy of Ooyama (1990). The equations are formulated using hydrostatic pressure as an
independent variable (Laprise, 1992). The vertical coordinate is terrain following, using a hybrid
σ − p formulation. In this chapter we define the vertical coordinate, present the moist flux-form
equations in Cartesian space, and augment the equations to include projections to the sphere.
η=1 ps
7
In the ARW Version, 4 the vertical coordinate has been generalized to allow the influence of
the terrain on the coordinate surfaces to be removed more rapidly with increasing height above
the surface, as illustrated in Fig. 2.1b. For this modified vertical coordinate, we employ a hybrid
sigma-pressure vertical coordinate as described by Park et al. (2013), which is similar to that
used in the National Center for Atmospheric Research (NCAR) Community Atmospheric Model
(CAM),as described in the NCAR CAM3.0 Technical Note:
pd = B(η)(ps − pt ) + [η − B(η)](p0 − pt ) + pt , (2.2)
where p0 is a reference sea-level pressure. (This coordinate representation differs somewhat
from CAM in that it is based on dry pressure instead of full pressure and is normalized using
pt such that η = 0 at pd = pt .) Here, B(η) defines the relative weighting between the terrain-
following sigma coordinate and a pure pressure coordinate, such that η corresponds to the sigma
coordinate (2.1) for B(η) = η and reverts to a hydrostatic pressure coordinate for B(η) = 0. To
smoothly transition from a sigma coordinate near the surface to a pressure coordinate at upper
levels, B(η) is defined by a third order polynomial
B(η) = c1 + c2 η + c3 η 2 + c4 η 3 (2.3)
(where the subscript η denotes differentiation) subject to the boundary conditions
B(1) = 1, Bη (1) = 1, B(ηc ) = 0, Bη (ηc ) = 0, (2.4)
such that
2ηc2 −ηc (4 + ηc + ηc2 ) 2(1 + ηc + ηc2 ) −(1 + ηc )
c1 = , c2 = , c3 = , c4 = , (2.5)
(1 − ηc )3 (1 − ηc )3 (1 − ηc )3 (1 − ηc )3
where ηc is a specified value of η at which it becomes a pure pressure coordinate.
Figure 2.2 displays the B(η) profiles for the traditional sigma coordinate and for the hybrid
coordinate for several values of the parameter ηc . Plotted as a function of η (Fig. 2.2a),
these profiles depict the form of the polynomial defined in (2.3). However, plotting B(η) as a
function of height (Fig. 2.2b) provides a better physical sense of the transition toward a pressure
coordinate with increasing height. For example, with a model domain having a depth of 30 km,
for ηc = 0.2 the vertical coordinate becomes a pure pressure coordinate at an altitude of about
12 km.
The vertical coordinate metric is defined as
∂pd
µd = = Bη (η)(ps − pt ) + [1 − Bη (η)](p0 − pt ). (2.6)
∂η
Since µd ∆η = ∆pd = −gρd ∆z is proportional to the mass per unit area within a grid cell, the
appropriate flux forms for the prognostic variables are defined as
V = µd v = (U, V, W ), Ω = µd ω, Θm = µd θm , Q m = µd q m . (2.7)
Here, v = (u, v, w) are the covariant velocities in the horizontal and vertical directions, while
ω = η̇ is the contravariant ‘vertical’ velocity. θm = θ(1 + (Rv /Rd )qv ) ≈ θ(1 + 1.61qv ) is the moist
potential temperature and Qm = µd qm , where qm = qv , qc , qr ... represents the mixing ratios of
moisture variables (water vapor, cloud water, rain water, ...). Although the geopotential φ = gz
is also a prognostic variable in the governing equations of the ARW, it is not written in flux-form
as µd φ is not a conserved quantity.
8
0 30
(a) (b)
σ
η
0 20
.2
z (km)
0.5 .4
σ
10 0
ηc .2
.4
ηc
1.0 0
0 0.5 1.0 0 0.5 1.0
B(η) B(η)
Figure 2.2: B(η) profiles for sigma (σ) coordinate and for hybrid coordinate for ηc =
0., 0.1, 0.2, 0.3, 0.4, and 0.5 displayed (a) as a function of η, and (b) as a function of height for
a standard atmosphere in a domain with a 30 km upper boundary.
9
and
V · ∇a = U ∂x a + V ∂y a + Ω∂η a,
where a represents a generic variable. γ = cp /cv = 1.4 is the ratio of the heat capacities for
dry air, Rd is the gas constant for dry air, and p0 is a reference surface pressure (typically 105
Pascals). The right-hand-side (RHS) terms FU , FV , FW , and FΘm represent forcing terms arising
from model physics, turbulent mixing, spherical projections, and the earth’s rotation.
In specifying the prognostic set of equations, a prognostic pressure equation could be used in
place of (2.13) (see Laprise, 1992). However, pressure is not a conserved variable and with pres-
sure as a prognostic variable we could not use the conservation equation (2.11) for Θm because
they are linearly dependent. Additionally, prognostic pressure equations have the disadvantage
of possessing a mass divergence term multiplied by a large coefficient (proportional to the sound
speed) which makes spatial and temporal discretization problematic. It should be noted that
the relation for the dry hydrostatic pressure (2.15) does not impose a hydrostatic constraint on
the solution, rather it is a diagnostic relation that formally is part of the coordinate definition.
In the hydrostatic counterpart to the nonhydrostatic equations, the full hydrostatic equation
∂η p = µd αd /α replaces the vertical momentum equation (2.10) and enforces a hydrostatic con-
straint on the solution.
In previous versions of the ARW the prognostic thermodynamic equation was expressed in
terms of Θ instead of Θm in (2.11). That representation is consistent with the expectation
that variations in qv have little impact during the small acoustic time steps in the split-explicit
numerics described in the next chapter, which has proven to be a robust approximation over a
wide spectrum of applications. However, in a recent study, Xiao et al. (2015) found that in
high-resolution LES with sharp vertical variations in water vapor, WRF simulations exhibited
a strong sensitivity to the time step used to accommodate the acoustic modes. Recasting the
prognostic thermodynamic equation in terms of Θm allowed consistent treatment of moisture
in the calculation of pressure during the acoustic time steps, and spurious motions and time-
step sensitivity were eliminated. Use of Θm as the prognostic thermodynamic variable is also
consistent with formulation used in the Model for Prediction Across Scales (MPAS), which
utilizes similar split-explicit numerics (Skamarock et al., 2012).
10
In the ARW’s computational space, ∆x and ∆y are constants. Orthogonal projections to
the sphere require that the physical distances between grid points in the projection vary with
position on the grid. To transform the governing equations, map scale factors mx and my are
defined as the ratio of the distance in computational space to the corresponding distance on the
earth’s surface:
(∆x, ∆y)
(mx , my ) = . (2.17)
distance on the earth
The ARW solver includes the map-scale factors in the governing equations by redefining the
momentum variables as
Using these redefined momentum variables, the governing prognostic equations (2.8)-(2.14) in-
cluding map factors can be written as
∂t U + mx [∂x (U u) + ∂y (V u)]
+ ∂η (Ωu) + (mx /my )[µd α∂x p + (α/αd )∂η p∂x φ] = FU (2.18)
∂t V + my [∂x (U v) + ∂y (V v)]
+(my /mx )∂η (Ωv) + (my /mx )[µd α∂y p + (α/αd )∂η p∂y φ] = FV (2.19)
∂t W + mx [∂x (U w) + ∂y (V w)] + ∂η (Ωw) − m−1y g[(α/αd )∂η p − µd ] = FW (2.20)
∂t Θm + mx my [∂x (U θm ) + ∂y (V θm )] + my ∂η (Ωθm ) = FΘm (2.21)
∂t µd + mx my [Ux + Vy ] + my ∂η (Ω) = 0 (2.22)
−1
∂t φ + µd [mx my (U ∂x φ + V ∂y φ) + my Ω∂η φ − my gW ] = 0 (2.23)
∂t Qm + mx my ∂x (U qm ) + ∂y (V qm )] + my ∂η (Ωqm ) = FQm , (2.24)
which are solved together with the diagnostic equations (2.15) and (2.16).
The right-hand-side terms of the momentum equations (2.18) – (2.20) contain the Coriolis
and curvature terms along with mixing terms and physical forcings. Including the map-scale
factors (2.17), the Coriolis and curvature terms are cast in the following form:
mx my ∂mx ∂my u
FUcor = + f +u −v V − + e cos αr W (2.25)
my mx ∂y ∂x re
my my ∂mx ∂my v
FVcor =− f +u −v U+ − e sin αr W (2.26)
mx mx ∂y ∂x re
mx 1 mx
FWcor = +e U cos αr − V sin αr + uU + vV , (2.27)
my re my
where αr is the local rotation angle between the y-axis and the meridians, ψ is the latitude,
f = 2Ωe sin ψ, e = 2Ωe cos ψ, Ωe is the angular rotation rate of the earth, and re is the radius of
the earth. In this formulation we have approximated the radial distance from the center of the
earth as the mean earth radius re , and we have not taken into account the change in horizontal
grid distance as a function of the radius. The terms containing the map-scale factors represent
11
the horizontal curvature terms, those containing re relate to vertical (earth-surface) curvature,
and those with e and f are the Coriolis force.
For the isotropic projections (Lambert conformal, polar stereographic, and Mercator), the
map-scale factors are the same in both horizontal directions, such that mx = my = m, where m
typically only varies with latitude. For the anisotropic latitude-longitude grid, mx = sec ψ and
my = 1, so that ∂mx /∂y = (1/re ) sec ψ tan ψ and ∂my /∂x = 0, given that ∂/∂y = (1/re )∂/∂ψ.
For idealized cases on a Cartesian grid, the map-scale factors mx = my = 1, f is specified, and
e and re−1 should be zero to remove the curvature terms.
and the mass conservation equation (2.22) and geopotential equation (2.23) become
and the diagnostic equation for dry hydrostatic pressure (2.15) becomes
Additionally, an option is available to use the hypsometric equation for the dry hydrostatic
pressure in place of (2.33):
∂φ/∂(ln pd ) = −pd (ᾱd + αd0 ). (2.34)
This form of the hydrostatic relation can produce a more accurate discretization compared with
(2.33) when the variation with η of temperature (pd αd ) is more linear than that of density
12
(αd−1 ). The conservation equations for the potential temperature (2.21) and the scalar moisture
equations (2.24) remain unchanged.
Equations (2.28) – (2.33), together with the equation of state (2.16), represent the equations
solved in the ARW. The RHS terms in these equations include the Coriolis terms (2.25) –
(2.27), mixing terms (described in Chapter 4), and parameterized physics (described in Chapter
8). Also note that the equation of state (2.16) cannot be written in perturbation form because
of the exponent in the expression. For small perturbation simulations, accuracy for perturbation
variables can be maintained by linearizing (2.16) for the perturbation variables.
13
14
Chapter 3
Model Discretization
∆t
Φ∗ = Φt + R(Φt ) (3.1)
3
∆t
Φ∗∗ = Φt + R(Φ∗ ) (3.2)
2
Φt+∆t = Φt + ∆tR(Φ∗∗ ) (3.3)
15
where ∆t is the time step for the low-frequency modes (the model time step). In (3.1) – (3.3),
superscripts denote time levels. This scheme is not a true Runge-Kutta scheme per se because,
while it is third-order accurate for linear equations, it is only second-order accurate for nonlinear
equations. With respect to the ARW equations, the time derivatives Φt are the partial time
derivatives (the leftmost terms) in equations (2.21), (2.24), and (2.28) – (2.32), and R(Φ) are
the remaining terms in those equations.
∗ ∗ ∗
V00 = V − Vt , Ω00 = Ω − Ωt , Θ00m = Θm − Θtm ,
t∗ t∗ t∗
φ00 = φ0 − φ0 , αd00 = αd0 − αd0 , µ00d = µ0d − µ0 d .
1 t∗ 00
αd00 00
= − t∗ ∂η φ + αd µd . (3.4)
µd
Additionally, we also introduce a version of the equation of state that is linearized about t∗ ,
∗ ∗
where c2s = γpt αdt is the square of the sound speed. The linearized state equation (3.5) and
the vertical coordinate definition (3.4) are used to cast the vertical pressure gradient in (2.30)
in terms of the model’s prognostic variables. By combining (3.5) and (3.4), the vertical pressure
gradient can be expressed as
c2s Θ00m
00 00
∂η p = ∂η (C∂η φ ) + ∂η ∗ , (3.6)
αdt Θtm∗
∗ ∗2
where C = c2s /µtd αdt . This linearization about the most recent large time step should be highly
accurate over the time interval of the several small time steps.
These variables along with (3.6) are substituted into equations (2.21) and (2.28) – (2.33) and
16
lead to the acoustic time-step equations:
∗ ∗ ∗ ∗ τ τ τ ∗ τ ∗
∂t U 00 + (mx /my )(αt /αdt ) µtd αdt ∂x p00 + αd00 ∂x p + ∂x φ00 + ∂x φt (∂η p00 − µ00d ) = RUt (3.7)
∗ ∗ ∗ ∗ τ τ τ ∗ τ ∗
∂t V 00 + (my /mx )(αt /αdt ) µtd αdt ∂y p00 + αd00 ∂y p + ∂y φ00 + ∂y φt (∂η p00 − µ00d ) = RVt (3.8)
t∗
δτ µ00d 00
+ mx my [∂x U + ∂y V ] 00 τ +∆τ
+ my ∂η Ω00τ +∆τ
= Rµ (3.9)
t∗ 00 t∗ τ +∆τ t∗ t∗
δτ Θ00m + mx my [∂x (U 00 θm
+ ∂y (V ) θm )] + my ∂η (Ω00τ +∆τ θm ) = RΘm (3.10)
2 00 τ
00 −1 cs Θm ∗
∗
δτ W − my g (α/αd )t ∂η (C∂η φ00 ) + ∂η − µ00d = RW t (3.11)
αt∗ Θtm∗
1 ∗ τ ∗
δτ φ00 + t∗
[my Ω00τ +∆τ δη φt − my gW 00 ] = Rφ t . (3.12)
µd
The RHS terms in (3.7) – (3.12) are fixed for the acoustic steps that comprise the time integration
of each RK3 sub-step (i.e., (3.1) – (3.3)), and are given by
∗
RUt = − mx [∂x (U u) + ∂y (V u)] − ∂η (Ωu)
− (mx /my )(α/αd ) [µd (∂x φ0 + αd ∂x p0 + αd0 ∂x p) + ∂x φ(∂η p0 − µ0d )] (3.13)
t∗
RV = − my [∂x (U v) + ∂y (V v)] − (my /mx )∂η (Ωv)
− (my /mx )(α/αd ) [µd (∂y φ0 + αd ∂y p0 + αd0 ∂y p) + ∂y φ(∂η p0 − µ0d )] (3.14)
∗
Rµt d = − mx my [∂x U + ∂y V ] − my ∂η Ω (3.15)
t∗
RΘ m
= − mx my [∂x (U θm ) + ∂y (V θm )] − my ∂η (Ωθm ) + FΘm (3.16)
t∗
RW = − mx [∂x (U w) + ∂y (V w)] − ∂η (Ωw)
+ m−1 0 −1 0
y g(α/αd )[∂η p − µ̄d (qv + qc + qr )] − my µd g + FW (3.17)
t∗
Rφ = − µ−1
d [mx my (U ∂x φ + V ∂y φ) + my Ω∂η φ − my gW ], (3.18)
where all variables in (3.13) – (3.18) are evaluated at time t∗ (i.e., using Φt , Φ∗ , or Φ∗∗ for the
appropriate RK3 sub-step in (3.1) – (3.3)). Equations (3.7) – (3.12) utilize the discrete acoustic
time-step operator
aτ +∆τ − aτ
δτ a = ,
∆τ
where ∆τ is the acoustic time step, and terms averaged in time over an acoustic time step are
slightly forward centered using an averaging operator
1 + β τ +∆τ 1 − β τ
aτ = a + a , (3.19)
2 2
where β is a user-specified parameter (see Section 4.3.3).
The integration over the acoustic time steps proceeds as follows. Beginning with the small
time-step variables at time τ , (3.7) and (3.8) are stepped forward to obtain U 00 τ +∆τ and V 00 τ +∆τ .
Both µ00 τ +∆τ and Ω00 τ +∆τ are then calculated from (3.9). This is accomplished by first integrating
(3.9) vertically from the surface to the material surface at the top of the domain, which removes
the ∂η Ω00 term. Recalling that µd = ∂pd /∂η and that pd for the hybrid vertical coordinate is
defined by (2.2), the vertical integral of (3.9) becomes,
Z 0
00
δτ pc = mx my [∂x U 00 + ∂y V 00 ]τ +∆τ dη, (3.20)
1
17
where pc (x, y) = ps − pt is the dry hydrostatic pressure difference (mass) in the vertical column
at (x, y). After computing δτ p00c from (3.20), Ω00 τ +∆τ is obtained by vertically integrating the
∂η Ω00 term (3.9) (with δτ µ = Bη δτ pc ) using the lower boundary condition Ω00 = 0 at the surface,
which yields Z η
00 τ +∆τ −1
[∂x U 00 + ∂y V 00 ]τ +∆τ dη,
Ω = my 1 − B(η) δτ pc − mx (3.21)
1
and µ00τ
d
+∆τ
is recovered using (2.2):
τ +∆τ τ +∆τ
µ00d (x, y, η) = Bη (η)p00c
(x, y) + 1 − Bη (η) (p0 − pt ). (3.22)
From (3.22), it is evident that µ need not be stored as a three-dimensional array, but can
be readily constructed when needed from the two-dimensional pc (x, y) array together with the
one-dimensional Bη (η) profile.
Knowing Ω00 τ +∆τ , (3.10) can be stepped forward to calculate Θ00m τ +∆τ . Equations (3.11) and
(3.12) are combined to form a vertically implicit equation that is solved for W 00 τ +∆τ subject to
the boundary condition Ω = Ω00 = 0 at the surface (z = h(x, y)) and p0 = 0 along the model
top. φ00 τ +∆τ is then obtained from (3.12), and p00 τ +∆τ and αd00 τ +∆τ are recovered from (3.5) and
(3.4).
18
Begin Time Step
Begin RK3 Loop: Steps 1, 2, and 3
(1) If RK3 step 1, compute and store FΦ
(i.e., physics tendencies for RK3 step, including mixing).
∗
(2) Compute Rt , (3.13)–(3.18)
Begin Acoustic Step Loop: Steps 1 → n,
RK3 step 1, n = 1, ∆τ = ∆t/3;
RK3 step 2, n = ns /2, ∆τ = ∆t/ns ;
RK3 step 3, n = ns , ∆τ = ∆t/ns .
(3) Advance horizontal momentum, (3.7) and (3.8)
Global: Apply polar filter to U 00τ +∆τ and V 00τ +∆τ .
(4) Advance µd (3.9) and compute Ω00 τ +∆τ then advance Θm (3.10)
Global: Apply polar filter to µτd+∆τ and Θm 00τ +∆τ
.
(5) Advance W and φ (3.11) and (3.12)
Global: Apply polar filter to W 00τ +∆τ and φ00τ +∆τ .
(6) Diagnose p00 and α00 using (3.5) and (3.4)
End Acoustic Step Loop
(7) Scalar transport: Advance scalars (2.24)
over RK3 substep (3.1), (3.2) or (3.3)
(using mass fluxes U , V and Ω time-averaged over the acoustic steps).
Global: Apply polar filter to scalars.
(8) Using updated prognostic variables, compute p0 with (2.16) and α0 with
(2.33) or (2.34)
End RK3 Loop
(9) Compute non-RK3 physics (currently microphysics), advance variables.
Global: Apply polar filter to updated variables.
End Time Step
Figure 3.1: Time step integration sequence. Here n represents the number of acoustic time
steps for a given substep of the RK3 integration, and ns is the ratio of the RK3 time step to
the acoustic time step for the second and third RK3 substeps.
19
3.1). Additional diabatic contributions are integrated in an additive-time-split manner in step
(9) after the RK3 update is complete. Thus, the diabatic forcing computed in step (9) (the
t∗
microphysics in the current release of the ARW) does not appear in RΘ m
from (3.10) in the
acoustic integration. We have discovered that this time splitting can excite acoustic waves and
can give rise to noise in the solutions for some applications. Note that the non-RK3 physics are
integrated in step (9) because balances produced in the physics are required at the end of the
time step (e.g., the saturation adjustment in the microphysics). So while moving these non-RK3
physics into step (1) would eliminate the noise, the balances produced by these physics would
be altered.
We have found that the excitation of the acoustic modes can be circumvented while leaving
the non-RK3 physics in step (9) by using the following procedure that is implemented in the
ARW. In step (1) of the integration procedure (Fig. 3.1), an estimate of the diabatic forcing in
the Θm equation arising from the non-RK3 physics in step (9) is included in the diabatic forcing
t∗
term RΘ m
in (3.10) (which is advanced in step 4). This estimated diabatic forcing is then
removed from the updated Θm after the RK3 integration is complete and before the evaluation
of the non-RK3 physics in step (9). We use the diabatic forcing from the previous time step
as the estimated forcing; hence this procedure results in few additional computations outside of
saving the diabatic forcing between time steps.
20
y η
vi,j+ 3/2 vi+1,j+ 3/2 wi,k+ 3/2 wi+1,k+ 3/2
{{
vi,j+ 1/2 vi+1,j+ 1/2
wi,k+ wi+1,k+
{
Δηk+1/2 1/2 1/2
vi,j-
1/2 vi+1,j- 1/2 wi,k- 1/2 wi+1,k- 1/2
x
{
{
Δx Δx
horizontal grid vertical grid
computed at mass points. The grid lengths ∆x and ∆y are constants in the model formulation;
changes in the physical grid lengths associated with the various projections to the sphere are
accounted for using the map factors introduced in Section 2.3. The vertical grid length ∆η is
not a fixed constant; it is specified in the initialization. The user is free to specify the η values
of the model levels subject to the constraint that η = 1 at the surface, η = 0 at the model
top, and η decreases monotonically between the surface and model top. Using these grid and
variable definitions, we can define the spatial discretization for the ARW solver.
µd u µd x u µd v µd y v
U= → , V = → ,
my my x mx mx y
where the discrete operator ax denotes a linear interpolation operator. The grid lengths ∆x and
∆y are constant, hence in this case the operator reduces to ax = (ai+1/2 + ai−1/2 )/2.
Using these definitions, we can write the spatially discrete acoustic step equations
21
(3.7) – (3.12) as
00 ∗ x t∗
x 00 τ
∗ x
00 τ x τη
∂t U + (mx /my )(αt∗ /αdt ) µtd
αd ∂x p + αd ∂x p + ∂x φ 00
τ
∗η xη 00 x ∗
+∂x φ t ∂η p 00 − µd = RUt (3.23)
y y y
00 t∗ t∗ t∗ t∗ 00 τ 00 τ y 00 τη
∂t V + (my /mx )(α /αd ) µd αd ∂y p + αd ∂y p + ∂y φ
τ
η y η y ∗
+∂y φt ∗
∂η p00 − µ00d = RVt (3.24)
∗
δτ µ00d + mx my [δx U 00 + δy V 00 ]τ +∆τ + my δη Ω00τ +∆τ = Rµt (3.25)
x 00 t∗ y τ +∆τ η t∗
δτ Θ00m + mx my [δx (U 00 θm
t∗ )
+ δy (V θm )] + my δη (Ω00τ +∆τ θmt∗ ) = RΘm (3.26)
2 00 τ
00 −1 ∗η cs Θm ∗
δτ W − my g (α/αd )t δη (Cδη φ00 ) + δη − µd00
= RW t (3.27)
αt∗ Θtm∗
1 η τ ∗
δτ φ00 + t∗
[my Ω00τ +∆τ δη φt∗ − my gW 00 ] = Rφ t , (3.28)
µd
where the discrete operator
δx a = ∆x−1 (ai+1/2 − ai−1/2 ) (3.29)
with the operators δy and δη similarily defined. Additionally, the operator aη is a vertical
interpolation operator. Using the notation given for the vertically stretched grid depicted in
Fig. 3.2, it is defined as
η 1 ∆ηk ∆ηk+1
a |k+1/2 = ak+1 + ak . (3.30)
2 ∆ηk+1/2 ∆ηk+1/2
This operator vertically interpolates variables on mass levels k to the w levels (k + 12 ). It should
be noted that the vertical grid is defined such that vertical interpolation from w levels to mass
levels reduces to aηk = (ak+1/2 + ak−1/2 )/2 (see Fig. 3.2).
The RHS terms in the discrete acoustic step equations for momentum (3.23), (3.24) and
(3.27) are discretized as
xη
x
h i
t∗ x 0
η x 0 0x η
0 0x
RU = −(mx /my )(α/αd ) µd (∂x φ + αd ∂x p + αd ∂x p) + ∂x φ (∂η p − µd )
+ FUcor + advection + mixing + physics, (3.31)
yη
y
h i
t∗ y 0
η y 0 0y η
0 0y
RV = −(my /mx )(α/αd ) µd (∂y φ + αd ∂y p + αd ∂y p) + ∂y φ (∂η p − µd )
+ FVcor + advection + mixing + physics, (3.32)
t∗ η
RW = m−1 0 η −1 0
y g(α/αd ) [δη p +µ̄d qm ] − my µd g
+ FWcor + advection + mixing + buoyancy + physics. (3.33)
22
map projections (Lambert conformal, polar stereographic, and Mercator) where mx = my = m,
their spatial discretization is
xη
x x y x uW
xy x xη x
FUcor = + f + u δy m − v δx m V − e W cos αr − , (3.34)
re
yη
y x y y xy y yη y vW
FVcor = − f + u δy m − v δx m U + e W sin αr − , (3.35)
re
xη xη yη
xη yη u U + v yη V
FWcor = +e(U cos αr − V sin αr ) + . (3.36)
re
xy xy xη yη
Here the operators () = () , and likewise for () and () .
For the non-isotropic latitude longitude projection, the Coriolis and curvature terms are
discretized as
xy xη
mx x uV xy xη uW
FUcor = f V + tan ψ − ex W cos αr x − , (3.37)
my re re
xy yη
uxy U
my y xy y yη y vW
FVcor = −f U − tan ψ + e W sin αr − , (3.38)
mx re re
xη xη yη
xη yη u U + (mx /my )v yη V
FWcor = +e(U cos αr − (mx /my )V sin αr ) + . (3.39)
re
3.2.3 Advection
The advection terms in the ARW solver are in the form of a flux divergence and are a subset of
the RHS terms in equations (3.13) – (3.18):
∗
RUt adv = − mx [∂x (U u) + ∂y (V u)] + ∂η (Ωu) (3.40)
∗
RVt adv = − my [∂x (U v) + ∂y (V v)] + (mx /my )∂η (Ωv) (3.41)
∗
Rµt adv = − mx my [∂x U + ∂y V ] + my ∂η Ω (3.42)
t∗
RΘ adv
= − mx my [∂x (U θm ) + ∂y (V θm )] − my ∂η (Ωθm ) (3.43)
t∗
RWadv = − mx [∂x (U w) + ∂y (V w)] + ∂η (Ωw) (3.44)
∗
Rφt adv =− µ−1
d [mx my (U ∂x φ + V ∂y φ) + my Ω∂η φ]. (3.45)
For the mass conservation equation, the flux divergence is discretized using a 2nd-order centered
approximation:
∗ ∗ ∗
Rµt adv = −mx my [δx U + δy V ]t + my δη Ωt . (3.46)
In the current version of the ARW, the advection of vector quantities (momentum) and scalars
is performed using the RK3 time integration as outlined in Fig. 3.1. The spatial discretization
used in this approach is outlined in the next section. For many applications it is desirable to
use positive definite or monotonic advection schemes for scalar transport.
23
RK3 Advection
2nd through 6th order accurate spatial discretizations of the flux divergence are available in the
ARW for momentum, scalars and geopotential using the RK3 time-integration scheme (scalar
advection option 1, step 7 in the time-split integration sequence in Fig. 3.1). The discrete
operators can be illustrated by considering the flux divergence equation for a scalar q in its
discrete form:
∗
Rqt adv = −mx my [δx (U q xadv ) + δy (V q yadv )] − my δη (Ωq ηadv ). (3.47)
As in the pressure gradient discretization, the discrete operator is defined as
δx (U q xadv ) = ∆x−1 (U q xadv )i+1/2 − (U q xadv )i−1/2 .
(3.48)
The different order advection schemes correspond to different definitions for the operator q xadv .
The even order operators (2nd , 4th , and 6th ) are
1
2nd order: (q xadv )i−1/2 = (qi + qi−1 )
2
th xadv 7 1
4 order: (q )i−1/2 = (qi + qi−1 ) − (qi+1 + qi−2 )
12 12
th xadv 37 2 1
6 order: (q )i−1/2 = (qi + qi−1 ) − (qi+1 + qi−2 ) + (qi+2 + qi−3 ),
60 15 60
rd th
and the odd order operators (3 and 5 ) are
th
3rd order: (q xadv )i−1/2 = (q xadv )4i−1/2
1
+ sign(U ) (qi+1 − qi−2 ) − 3(qi − qi−1 )
12
th
5th order: (q xadv )i−1/2 = (q xadv )6i−1/2
1
− sign(U )(qi+2 − qi−3 ) − 5(qi+1 − qi−2 ) + 10(qi − qi−1 ) .
60
The even-order advection operators are spatially centered and thus contain no implicit dif-
fusion outside of the diffusion inherent in the RK3 time integration. The odd-order schemes are
upwind-biased, and the spatial discretization is inherently diffusive. The behavior of the up-
wind schemes is easily understood by expanding (3.48) using the 5th order operator, assuming
a constant mass flux U and multiplying by the timestep ∆t:
U ∆t 1
∆tδx (U q xadv ) = ∆tδ(U q)|6th − (−qi−3 + 6qi−2 − 15qi−1 + 20qi − 15qi+1 + 6qi+2 − qi+3 )
∆x 60
Cr 6 ∂ 6 q
= ∆tδ(U q)|6th − ∆x + higher order terms
60 ∂x6
Similarly, we can expand (3.48) using the 3rd order operator:
Cr 4 ∂ 4 q
∆tδx (U q xadv ) = ∆tδ(U q)|4th + ∆x + higher order terms
12 ∂x4
As is evident in their formulation, the odd-order schemes are comprised of the next higher (even)
order centered scheme plus an upwind term that, for a constant transport mass flux, is a diffusion
term of that next higher (even) order with a hyper-viscosity proportional to the Courant number
(Cr). The diffusion term is the leading order error term in the flux divergence discretization.
Further details concerning RK3 advection can be found in Wicker and Skamarock (2002)
24
Positive-Definite and Shape-Preserving Limiters for RK3 Advection
Mixing ratios of moisture, chemical species or other tracer species should remain positive-
definite, that is, negative masses should not be permitted. Additionally, transport should
not introduce any new maxima or minima in the mixing ratios; this property is referred to
as shape-preservation also monotonicity. The Runge-Kutta transport integration defined by
the timestepping algorithm (3.1) – (3.3), combined with the flux divergence operator (3.47), is
conservative but it does not guarantee positive definiteness or shape preservation; any negative
values will be offset by positive mass such that mass is conserved, and new maxima and minima
may be produced. In many physics options, negative mixing ratios will be set to zero, and this
will result in an increase in mass of that species. Both positive-definite and shape-preserving
flux renormalizations, applied on the final Runge-Kutta transport step (3.1), can be used to
remove these unphysical effects from the RK3 scalar transport scheme. The positive-definite
renormalization in the ARW solver is described in Skamarock and Weisman (2008) and the
shape-preserving (monotonic) extension is described in Wang et al. (2009), both of which follow
an approach first described by Zalesak (1979).
The renormalizations occur as part of the final RK3 time integration step for scalar transport,
and the preliminary (i.e. before renormalization) evaluation of the final transport step can be
expressed as
where the flux divergence is evaluated using the (**) time level predicted in RK3 step (3.2).
The positive-definite flux renormalization replaces (3.2) with the following two steps. First, the
scalar mixing ratio is updated using the tendency derived from the model physics and source/sink
terms.
(µφ)∗∗∗ = (µφ)t + ∆tµSφt (3.50)
where we denote this new predictor as (µφ)∗∗∗ . Second, the full update is computed using a flux
divergence composed of a first-order upwind flux plus a higher order correction:
In (3.51), ()1 denotes a first-order upwind flux and R()0 denotes a renormalized higher-order
correction flux. The higher-order correction flux is the difference between the full RK3 flux and
the first-order upwind flux, that is,
with similar definitions for (V q yadv )0 and (Ωq ηadv )0 . The correction flux can then be renormalized
to achieve positive-definite or shape preserving (monotonic) behavior as follows. First, the
upwind fluxes are used to perform a partial update of the scalar mass.
˜ = (µφ)∗∗∗ − ∆t mx my [δx (U q xadv )1 + δy (V q yadv )1 ] + my δη (Ωq ηadv )1
(µφ)
25
This update is positive definite (and monotonic) by design because it is a property of the the
first-order upwind scheme. Next, a prediction of the minimum and maximum possible values of
the new-time-level species mass is computed at each point by using only the outward directed
fluxes (fluxes that remove mass from the control volume),
(µφ)t+∆t ˜ xadv 0
) + δy (V+ q yadv )0 ] + my δη (Ω+ q ηadv )0
min = (µφ) − ∆t mx my [δx (U+ q (3.53)
where U+ , V+ and Ω+ indicated the use of fluxes out of the control volume only, that is, only
those that contribute to lowering the scalar mass, and U− , V− and Ω− indicated the use of fluxes
into the control volume only, that is, only those that contribute to rasing the scalar mass.
For the positive-definite renormalization, from (3.53) the scalar mass (µφ)t+∆t min < 0 if and
only if
˜ < ∆t mx my [δx (U+ q xadv )0 + δy (V+ q yadv )0 ] + my δη (Ω+ q ηadv )0 .
(µφ)
For each volume where negative mass is indicated by (3.53), the fluxes are renormalized such
that the outgoing fluxes and mass in the volume are equivalent.
˜
µφ
R (U+ q xadv )0 = (U+ q xadv )0
∆t mx my [δx (U+ q xadv )0 + δy (V+ q yadv )0 ] + my δη (Ω+ q ηadv )0
with a similar renormalization applied to the (V+ q yadv )0 and (Ω+ q ηadv )0 .
For the shape-preserving (monotonic) renormalization, the minimum (maximum) value for a
cell mixing ratio at the new time level is set as the minimum (maximum) mixing ratio at time t
from that cell and its nearest neighbors. This lower bound on the mixing ratio, denoted (φ)min ,
will be violated when
˜ − (µφ)min
µφ
R (U+ q xadv )0 = (U+ q xadv )0
∆t mx my [δx (U+ q xadv )0 + δy (V+ q yadv )0 ] + my δη (Ω+ q ηadv )0
with a similar renormalization applied to the (V+ q yadv )0 and (Ω+ q ηadv )0 . A similar limiting process
is used to renormalize the fluxes such that the maximum mixing ratio is not exceeded. Note
that outgoing fluxes for one cell are inflow fluxes for another, and the renormalized flux used in
the update is the minimum flux necessary to satisfy the mimimum or maximum bounds for the
two cells utilizing the flux.
The renormalized higher-order-correction fluxes along with the first order fluxes are then
used in the update equation (3.51). Note that if no renormalization is needed the scheme (3.51)
reverts to the standard RK3 update because of the definition of the higher-order correction
(3.52).
26
U
V
µ, W
re∆ψ
re∆ψ
re∆ψ
Figure 3.3: Latitude-longitude grid structure in the pole region. In the ARW formulation,
re ∆ψ = ∆y/my .
27
in the acoustic sub-steps of the time-split integration procedure, see Section 3.1.2). Both are
limited by Courant numbers. In the following sections we describe how to choose time steps for
applications.
Table 3.1: Maximum stable Courant numbers for one-dimensional linear advection. From Wicker
and Skamarock (2002).
As is indicated in the table, the maximum stable Courant numbers for advection in the RK3
scheme are almost a factor of two greater than those for the leapfrog time-integration√scheme.
For advection in three spatial dimensions, the maximum stable Courant number is 1/ 3 times
the Courant numbers given in Table 3.1. For stability, the time step used in the ARW should
produce a maximum Courant number less than that given by theory. Thus, for 3D applications,
the time step should satisfy the following equation:
Crtheory ∆x
∆tmax < √ · , (3.55)
3 umax
where Crtheory is the Courant number taken from the RK3 entry in Table 3.1 and umax is the
maximum velocity expected in the simulation. For example in real-data applications, where jet
stream winds may reach as high as 100 ms−1 , the maximum time step would be approximately
80 s on a ∆x = 10 km grid using 5th order advection. For convection-permitting resolutions
(typically ∆x ≤ 5 km), the vertical velocities in convective updrafts produce the stability-
limiting Courant numbers. Given the additional constraint from the time splitting, and to
provide a safety buffer, we usually choose a time step that is approximately 25% less than that
given by (3.55). This time step is typically a factor of two greater than that used in leapfrog-
based models. For those users familiar with the MM5 model, the rule of thumb for choosing a
time step is that the time step, in seconds, should be approximately 3 times the horizontal grid
distance, in kilometers. For the ARW, the time step (in seconds) should be approximately 6
times the grid distance (in kilometers).
28
speed of √
sound. We typically use a more conservative estimate for this by replacing the limiting
value 1/ 2 with 1/2. Thus, the acoustic time step used in the model is
1 ∆x
∆τ < · . (3.56)
2 cs
For example, on a ∆x = 10 km grid, using a sound speed cs = 300 ms−1 , the acoustic time
step given in (3.56) is approximately 17 s. In the ARW, the ratio of the RK3 time step to the
acoustic time step must be an even integer. For our example using a ∆x = 10 km grid in a
real-data simulation, we would specify the RK3 time step ∆t = 60s (i.e., 25% less than the 80 s
step given by (3.55), and an acoustic time step ∆τ = 15 s (i.e., 1/4 of the RK3 step, rounding
down the ∆τ = 17 s step given by (3.56)). Note that it is the ratio of the RK3 time step to the
acoustic time step that is the required input in the ARW.
29
∆t = min(max(∆t, ∆tmin ), ∆tmax ), (3.62)
to guard the time step from becoming too small or too large. The computation of the number
of acoustic time steps ns , described in section 3.3.2, is handled during the first of each of the
RK3 model integration steps:
∆t
ns = max 2 · b300 · + 1c, 4 . (3.63)
∆x
For a simulation using nests, the fine-grid domain must maintain an integer number of
time steps within a single model integration step from the parent. At each of the starting model
integration steps for the child grid, those integration steps when the time on the fine-grid domain
equals the time on the parent’s domain, the adaptive time step algorithm is conducted for the
fine grid. Note that the ratio of the nominal grid distance between the parent and the child
does not necesarily imply that same ratio between the model integration time steps.
30
Chapter 4
A number of formulations for turbulent mixing and filtering are available in the ARW solver.
Some of these filters are used for numerical reasons. For example, divergence damping is used to
filter acoustic modes from the solution and polar filtering is used to reduce the timestep restric-
tion arising from the converging gridlines of the latitude-longitude grid. Other filters are meant
to represent sub-grid turbulence processes that cannot be resolved on the chosen grid. These
filters typically remove energy from the solution and are formulated in part on turbulence theory
and observations, or represent energy sink terms in some approximation to the Euler equation.
In this section, we begin by outlining the formulation and discretization of turbulent mixing
processes in the ARW solver commonly associated with sub-gridscale turbulence as parameter-
ized in cloud-scale models— the second-order horizontal and vertical mixing. In existing global
models and most NWP models, vertical mixing is parameterized within the planetary boundary
layer (PBL) physics. Vertical mixing parameterized within the PBL physics is described later in
Chapter 8. Here we note that, when a PBL parameterization is used, all other vertical mixing
is disabled. Following the outline of turbulent mixing parameterizations in this chapter, other
numerical filters available in the ARW solver are described.
where φ̂(k) and φ̂(k)f iltered are the Fourier coefficients for the generic variable φ before and after
31
filtering, and a(k) are the filter coefficients defined as a function of the dimensionless wavenumber
k. The ARW polar filter coefficients a(k) for the Fourier amplitudes are
" 2 !#
cos ψ 1
a(k) = min 1., max 0., 2
cos ψo sin (πk/n)
where ψ is the latitude on the computational grid, ψo is the latitude above which the polar filter
is applied (the filter is applied for |ψ| > ψo ), and n is the number of grid points in the latitude
circle. For even n, the grid admits wavenumber k = 0 and dimensionless wavenumbers ±k for
k = 1 → n/2 − 1, and the 2∆ wave k = n/2. For odd n, the grid admits wavenumber k = 0,
and dimensionless wavenumbers ±k for k = 1 → (n − 1)/2.
Polar filter applications within the split-explicit RK3 time integration scheme are given in
Figure 3.1. Within the acoustic steps, the filter is applied to U 00 τ +∆τ and V 00 τ +∆τ immediately
after the horizontal momentum is advanced, to pcτ +∆τ and Θ00m τ +∆τ immediately after the column
mass and the potential temperature are advanced, and to W 00 τ +∆τ and φ00 τ +∆τ after the vertically
implicit part of the acoustic step is complete. The polar filter is also applied to the scalars after
they are advanced every RK3 sub-step (3.1)-(3.3) and after the microphysics step.
The prognostic variables are coupled with the column mass µ when filtered, and in this way
mass is conserved. An exception to this is the geopotential φ which is filtered without being
coupled. The Fourier filtering is conservative, but it is not monotonic or positive definite. As
noted in the ARW stability discussion in Section 3.3.4, timestep stability should be based on the
longitudinal gridlength where the polar filters are activated, (∆x/mx )|ψo . The positive definite
transport option, presented in Section 3.2.3, will not necessarily produce positive-definite results
because of the converging gridlines and should not be used.
32
For the horizontal and vertical momentum equations, (4.1) is spatially discretized as
xy xη
∂t U = ... + µd x mx x δx (mx Kh δx u) + δy (my xy K h δy u) + m−1 2 x x −1 xη −1
y g (µd α ) δη (K v (α ) δη u),
xy yη
∂t V = ... + µd y my y δx (mx xy K h δx v) + δy (my Kh δy v) + m−1 2 y y −1 yη −1
x g (µd α ) δη (K v (α ) δη v),
xη yη
∂t W = ... + µd mx δx (mx x K h δx w) + δy (my y K h δy w) + m−1 2 η −1 −1
y g (µd α ) δη (Kv α δη w).
In the current ARW formulation for mixing on coordinate surfaces, the horizontal eddy viscosity
Kh is allowed to vary in space, whereas the vertical eddy viscosity does not vary in space; hence
there is no need for any spatial averaging of Kv . Additionally, note that the horizontal eddy
viscosity Kh is multiplied by the inverse turbulent Prandtl number Pr−1 for horizontal scalar
mixing.
Continuous Equations
The continuous equations for evaluating diffusion in physical space, using the velocity stress
tensor, are as follows for horizontal and vertical momentum:
∂t U = ... − mx ∂x τ11 + ∂y τ12 − zx ∂z τ11 − zy ∂z τ12 − ∂z τ13 (4.2)
∂t V = ... − my ∂x τ12 + ∂y τ22 − zx ∂z τ12 − zy ∂z τ22 − ∂z τ23 (4.3)
∂t W = ... − my ∂x τ13 + ∂y τ23 − zx ∂z τ13 − zy ∂z τ23 − ∂z τ33 . (4.4)
The stress tensor τ can be written as follows:
τ11 = −µd Kh D11
τ22 = −µd Kh D22
τ33 = −µd Kv D33
τ12 = −µd Kh D12
τ13 = −µd Kv D13
τ23 = −µd Kv D23 .
33
Symmetry sets the remaining tensor values; τ21 = τ12 , τ31 = τ13 , and τ32 = τ23 . The stress tensor
τ is calculated from the deformation tensor D. The continuous deformation tensor is defined as
D11 = 2 mx my ∂x (m−1 −1
y u) − zx ∂z (my u)
D22 = 2 mx my ∂y (m−1 −1
x v) − zy ∂z (mx v)
D33 = 2 ∂z w
D12 = mx my ∂y (m−1 −1 −1 −1
y u) − zy ∂z (my u) + ∂x (mx v) − zx ∂z (mx v)
D13 = mx my ∂x (m−1 −1
y w) − zx ∂z (my w) +∂z (u)
D23 = mx my ∂y (m−1 −1
y w) − zy ∂z (my w) +∂z (v).
The deformation tensor is symmetric, hence D21 = D12 , D31 = D13 , and D32 = D23 .
The diffusion formulation for scalars is
∂t (µd q) = ... + mx ∂x − ∂z zx µd mx Kh (∂x − zx ∂z ) +
my ∂y − ∂z zy µd my Kh (∂y − zy ∂z ) + ∂z µd Kv ∂z q. (4.5)
Spatial Discretization
Using the definition of the stress tensor, the spatial discretization of the ARW physical-space
diffusion operators for the horizontal and vertical momentum equations (4.2) - (4.4) are
∂t U = ... − mx x δx τ11 + δy τ12 − zx η δz τ11 xη − zy xyη δz τ12 yη − δz τ13
The discrete forms of the stress tensor and deformation tensor are
τ11 = −µd Kh D11
τ22 = −µd Kh D22
τ33 = −µd Kv D33
xy
τ12 = −µd xy Kh D12
xη
τ13 = −µd x Kv D13
yη
τ23 = −µd y Kv D23 ,
and
x x
u) − zx xη δz (m−1
D11 = 2 mx my δx (m−1
y y u)
y xη y
D22 = 2 mx my δy (m−1
x v) − zy δz (mx v)
−1
D33 = 2 δz w
xy x xη x yη y yη y xη
−1 −1 −1 −1
D12 = (mx my ) δy (my u) − zy δz (my u) + δx (mx v) − zx δz (mx v)
−1 xη
D13 = mx my δx (my w) − zx δz (m−1 w) +δz u
−1 yη
D23 = mx my δy (my w) − zy δz (m−1 y w) +δz v.
34
The spatial discretization for the scalar diffusion (4.5) is
xη
∂t (µd q) = ... + mx δx µd x H1 (q) − µd zx xη δz H1 (q)
yη
+ my δy µd y H2 (q) − µd zy yη δz H2 (q)
η
+ µ d δz K v δz q ,
where
x
H1 (q) = mx x Kh δx q − zx δz (q xη ) ,
y
H2 (q) = my y Kh δy q − zy δz (q yη ) .
The deformation tensor components have been defined in the previous section. The length scale
l = (∆x∆y)1/2 and Cs is a constant with a typical value Cs = 0.25. For scalar mixing, the
eddy viscosity is divided by the turbulent Prandtl number Pr that typically has a value of 1/3
(Deardorff, 1972). This option is most often used with a planetary boundary layer scheme that
independently handles the vertical mixing.
3D Smagorinsky Closure
The horizontal and vertical eddy viscosities can be determined using a 3D Smagorinsky turbu-
lence closure. This closure specifies the eddy viscosities as
2 2 2
1
−1 2 2
Kh,v = Cs lh,v max 0., D − Pr N , (4.6)
where
21 2 2 2 xy 2 xη 2 yη 2
D = D11 + D22 + D33 + D12 + D13 + D23 ,
2
and N is the Brunt-Väisälä frequency; the computation of N , including moisture effects, is
outlined in Section 4.2.4.
35
Two options are available for calculating the mixing length lh,v in (4.6). An isotropic length
scale (appropriate for ∆x, ∆y ' ∆z) can be chosen where lh,v = (∆x∆y∆z)1/3 and thus Kh =
Kv = K. The √ anisotropic option (appropriate for ∆x, ∆y >> ∆z) sets the horizontal mixing
length lh = ∆x∆y in the calculation of the horizontal eddy viscosity Kh using (4.6), and
lv = ∆z for the calculation of the vertical eddy viscosity Kv using (4.6).
Additionally, the eddy viscosities for scalar mixing are divided by the turbulent Prandtl
number Pr = 1/3.
Shear Production
The shear production term in (4.7) can be written as
2 2 2 xy xη yη
2 2 2
shear production = Kh D11 + Kh D22 + Kv D33 + Kh D12 + Kv D13 + Kv D23 .
36
Buoyancy
The buoyancy term in the TKE equation (4.7) is written as
buoyancy = −Kv N 2 ,
where the Brunt-Väisälä frequency N is computed using either the formula for a moist saturated
or unsaturated environment:
N 2 = g A∂z θe − ∂z qw
if qv ≥ qvs or qc ≥ 0.01 g/Kg;
2 1
N = g ∂z θ + 1.61∂z qv − ∂z qw if qv < qvs or qc < 0.01 g/Kg.
θ
where qw represents the total water (vapor + all liquid species + all ice species), L is the latent
heat of condensation and is the ratio of the molecular weight of water vapor to the molecular
weight of dry air. θe is the equivalent potential temperature and is defined as
Lqvs
θe = θ 1 + ,
Cp T
Dissipation
The dissipation term in (4.7) is
Ce3/2
dissipation = − ,
l
where
max(0, 0.93 − 1.9 Ck ) l
C = 1.9 Ck + ,
∆s
∆s = (∆x∆y∆z)1/3 ,
and
√
l = min (∆x∆y∆z)1/3 , 0.76 e/N .
37
slope was introduced. For the staggered velocities u and v, and for the cell centered variables
a, the filter can be expressed as
2−6 h 6 i
∆x mx x δx µd Fx (u) + ∆y 6 my x δy µd xy Fy (u) ,
∂t (µd u) = · · · β (4.8)
2∆t
2−6 h 6 i
∆x mx y δx µd xy Fx (v) + ∆y 6 my y δy µd Fy (v) ,
∂t (µd v) = · · · β (4.9)
2∆t
2−6 h 6 x
6 y
i
∂t (µd a) = · · · β ∆x mx δx µd Fx (a) + ∆y my δy µd Fy (a) . (4.10)
2∆t
β is a user-specified filtering coefficient representing the damping applied to 2∆ waves for each
filter application. The diffusive fluxes ∆x5 Fx = ∆x5 δx5 (a), i.e.
n o
∆x5 Fx (a)
= γs 10 a(i + 1) − a(i) − 5 a(i + 2) − a(i − 1) + a(i + 3) − a(i − 2) .
i+1/2
γs is a flux limiting factor that reduces the diffusive flux where the terrain slope is steep. γs = 1
if the limiter is inactive. When activated by the user, the flux limiting factor
−1 !
∂z ∂z
γs = max 1 − ,0 .
∂x ∂x cr
(dz/dx) is the slope of the coordinate surface, and it is approximated using the reference state
geopotential in WRF. (dz/dx)cr is a user specified critical value of the coordinate surface slope
at which the diffusive flux is completely removed. By design, 0 ≤ γs ≤ 1.
The user can choose to enforce monotonicity in the filtering. In the monotonic option, the
diffusive fluxes are set to zero if the diffusion is up-gradient, e.g. Fx (a) = 0 if Fx (a) δx (a) < 0
and Fy (a) = 0 if Fy (a) δy (a) < 0.
Note that the diffusive fluxes are computed in computational space; map factors are not
taken into account except in the flux divergence terms in (4.8) - (4.10) where they serve to
guarantee scalar mass conservation. The filter parameter β is dimensionless, thus the filter
scales with the timestep and gridsize. This explicit diffusion acts on all three components of
wind, potential temperature, all moisture variables and passive scalars, and subgrid turbulence
kinetic energy. Further details can be found in Knievel et al. (2007).
38
accomplished by forward weighting the pressure update in the acoustic step loop described in
Section 3.1.3, step (6). The linearized equation of state (3.5) is used to diagnose the pressure at
the new time τ after the U 00 , V 00 , µ00d , and Θ00m have been advanced. Divergence damping consists
of using a modified pressure in the computation of the pressure gradient terms in the horizontal
momentum equations in the acoustic steps (Equations (3.7) and (3.8)). Denoting the updated
value as pτ , the modified pressure p∗τ used in (3.7) and (3.8) can be written as
p∗τ = pτ + γd pτ − pτ −∆τ ,
(4.11)
∆x2
δτ U 00 = ... − γe δx (δτ −∆τ µ00d ) (4.12)
∆τ
and
∆y 2
δτ V 00 = ... − γe δy (δτ −∆τ µ00d ). (4.13)
∆τ
The quantity δτ −∆τ µ is the vertically-integrated mass divergence (i.e., (3.20)) from the previous
acoustic step (that is computed using the time τ values of U and V ), and γe is the external
mode damping coefficient. An external mode damping coefficient γe = 0.01 is typically used in
the ARW applications, independent of time step or grid size.
39
4.4.1 Absorbing Layer Using Spatial Filtering
This formulation of the absorbing layer increases the second-order horizontal and vertical eddy
viscosities in the absorbing layer using the following formulation:
∆x2
π ztop − z
Kdh = γg cos ,
∆t 2 zd
and
∆z 2
π ztop − z
Kdv = γg cos .
∆t 2 zd
Here γg is a user-specified damping coefficient, ztop is the height of the model top for a particular
grid column, zd is the depth of the damping layer (from the model top), and Kdh and Kdv are
the horizontal and vertical eddy viscosities for the wave absorbing layer. If other spatial filters
are being used, then the eddy viscosities that are used for the second-order horizontal and
vertical eddy viscosities are the maximum of (Kh , Kdh ) and (Kv , Kdv ). The effect of this filter
on gravity waves is discussed in Klemp and Lilly (1978), where guidance on the choice of the
damping coefficient γg can also be found.
where W̃ 00τ +∆τ is the solution to (3.11). The geopotential is then updated as usual using (3.12)
with the updated damped vertical velocity from (4.14). The perturbation pressure and specific
volume are computed as before using (3.5) and (3.4).
The variable τ (z) defines the vertical structure of the damping layer, and has a form similar
to the Rayleigh damper in Durran and Klemp (1983) and also used in the traditional Rayleigh
damping formulation discussed in Section 4.4.3
( h i
γr sin2 π2 1 − ztopzd−z for z ≥ (ztop − zd );
τ (z) =
0 otherwise,
40
where γr is a user specified damping coeficient, ztop is the height of the model top for a particular
grid column, and zd is the depth of the damping layer (from the model top). A typical value
for the damping coefficient for this formulation is γr = 0.2 s−1 (∼ 10N for typical stratospheric
values of the Brunt-Väisällä frequency). A complete analysis of this filter can be found in Klemp
et al. (2008).
∂t u = −τ (z) (u − u) ,
∂t v = −τ (z) (v − v) ,
∂t w = −τ (z)w,
∂t θ = −τ (z) θ − θ .
Overbars indicate the reference state fields, which are functions of z only and are typically
defined as the initial horizontally homogeneous fields in idealized simulations. The reference
state vertical velocity is assumed to be zero. The vertical structure of the damping layer is the
same as that used for the implicit vertical velocity damping described in Section 4.4.2. Because
the model surfaces change height with time in the ARW solver, the reference state values at
each grid point need to be recalculated at every time step. Thus, a linear interpolation scheme
is used to calculate updated reference state values based on the height of the model levels at
each time step.
The effect of this filter on gravity waves is discussed in Klemp and Lilly (1978), where
guidance on the choice of the damping coefficeint γr can also be found.
41
If Cr > Crβ , then
∂t W = ... − µd sign(W )γw (Cr − Crβ ),
where γw is the damping coefficient (typically 0.3 ms−2 ), and Crβ is the activation Courant
number (typically 1.0). The ARW outputs the location of this damping when it is active.
42
Chapter 5
Initial Conditions
ARW may be run with user-defined initial conditions for idealized simulations, or it may be run
using interpolated data from either an external analysis or forecast for real-data cases. Both 2D
and 3D test cases are provided for idealized simulations. Several sample cases are provided for
real-data simulations, which rely on pre-processing from an external package (usually the WRF
Preprocessing System, referred to as WPS) that converts the large-scale GRIB-formatted data
into an intermediate format suitable for ingest by the ARW’s real-data processor.
The programs that generate the specific initial conditions for the selected idealized or real-
data case function similarly. They provide ARW with:
• input data that are correctly staggered for the horizontal and vertical grid;
• hydrostatically balanced reference state and perturbation fields; and
• metadata specifying such information as the date, grid physical characteristics, and pro-
jection details.
Neither idealized nor real-data cases use observations to enhance their initial conditions; how-
ever, output from the ARW system initial condition program is suitable as input to the WRF
Variational Data Assimilation package (see Chapter 11).
43
Table 5.1: Ideal Cases. Listed are the available idealized cases for the Advanced Research WRF.
1D 2D 3D
Initial conditions for the Held-Suarez test case are analytically defined in the initialization
routines. Other than the baroclinic wave case that utilizes a binary file providing a 2D profile
specified in [y, z], the remaining idealized ARW cases use as input a 1D sounding specified
as a function of geometric height z. The text file defining this 1D sounding may be directly
edited by the user. The 1D specification of the sounding in these test files requires surface
values of pressure, potential temperature, and water vapor mixing ratio, followed by potential
temperature, vapor mixing ratio, and horizontal wind components at some heights above the
surface. Initialization programs for each case assume that this moist sounding represents an
atmosphere in hydrostatic balance.
Two sets of thermodynamic fields are needed for the model— the reference state and the
perturbation state (see Chapter 2 for discussion of the equations). The reference state used in
idealized initializations is computed using the input sounding from which moisture is discarded
(because the reference state is dry). The perturbation state is computed using the full moist
input sounding. The procedure for computing the hydrostatically-balanced ARW reference state
and perturbation state variables from the input sounding are as follows. First, density and both
a dry and full hydrostatic pressure are computed from the input sounding at the input sounding
levels z. This is accomplished by integrating the hydrostatic equation vertically up the column
using the surface pressure and potential temperature as the lower boundary condition. The
hydrostatic equation is
δz p = −ρd z g(1 + q v z ), (5.1)
where ρd z is a two-point average between input sounding levels, and δz p is the difference of
the pressure between input sounding levels, divided by the height difference. Additionally, the
equation of state is needed to close the system:
− ccv
1 Rd θm p p
αd = = , (5.2)
ρd po po
where qv and θm are given in the input sounding. (5.1) and (5.2) are a coupled set of nonlinear
equations for p and ρ in the vertical integration, solved by iteration. Dry pressure on input
sounding levels is computed by integrating the hydrostatic relation down from the top, excluding
the vapor component.
44
While the full pressure (dry, plus vapor) and dry air pressure are computed on the input
sounding levels, the pressure at the model top pt is computed by linear interpolation in height
(or possibly extrapolation), using the height of the model top (which is an input variable). The
column mass pc = ps − pt is computed by interpolating the dry air pressure to the surface, and
then subtracting pt . Given the column mass, the dry-air pressure at each η level is known from
the hybrid coordinate definition (2.2), and the coordinate metric µd is computed from (2.6).
Potential temperature from the input sounding is interpolated to each of the model pressure
levels, and the equation of state (5.2) is used to compute the inverse density αd . Finally, ARW’s
dry hydrostatic relation (2.15) is used to compute the geopotential. This procedure is used to
compute the reference state (based on a dry atmosphere) and the full state (using the full moist
sounding). Perturbation variables are computed as the difference between the full state and the
reference state. It should also be noted that in the nonhydrostatic model integration, inverse
density αd is diagnosed from the geopotential using this dry hydrostatic relation, and pressure
is diagnosed from the equation of state using the inverse density αd and the prognostic potential
temperature θm . Thus, ARW’s prognostic variables µd , θm , and φ are in exact hydrostatic
balance (to machine roundoff) with the model equations.
45
STATIC GEOGRID
DATA
REAL DATA
ARW
GRIB SYSTEM
UNGRIB REAL
DATA
Figure 5.1: Schematic showing data flow and program components in WPS, and how WPS feeds
initial data to ARW. Text in the rectangular boxes indicates program names. GEOGRID: defines
the model domain and creates static files of terrestrial data. UNGRIB: decodes GRIB-formatted
data. METGRID: interpolates meteorological data to the model domain.
type, land/water mask, map scale factors, map rotation angle, soil texture category, vegetation
greenness fraction, annual mean temperature, and latitude/longitude. After WPS processing,
the 2-dimensional time-dependent fields from the external model include surface pressure and
sea-level pressure (Pa), layers of soil temperature (K) and soil moisture (kg/kg, either total
moisture, or binned into total and liquid content), snow depth (m), skin temperature (K), sea
surface temperature (K), and a sea ice flag.
46
From (5.3), the 3-dimensional reference pressure (dry hydrostatic pressure pd ) is computed as a
function of the vertical coordinate η levels and the model top pt :
With (5.4), the reference temperature, though not permitted to get colder than Tmin , is a straight
line on a skew-T plot, defined as
pd
T d = max Tmin , T0 + A ln .
p0
For vertical locations where the pd < pstrat , the reference profile warms:
pd
T d = Tmin + γstrat ln .
pstrat
The isobaric temperature and the stratospheric correction supply a reasonable reference tem-
perature up to approximately 100 Pa. From the reference pressure and the final reference
temperature, the reference potential temperature is then defined as
RCd
p0 p
θd = T d , (5.5)
pd
and the reciprocal of the reference density using (5.4) and (5.5) is given by
− CCv
1 Rd θd pd p
αd = = . (5.6)
ρd p0 p0
∂ p̄d
µ̄d = = Bη (η)(p̄s − pt ) + 1 − Bη (η) (p0 − pt ). (5.7)
∂η
From (5.6) and (5.7), the reference state geopotential defined from the hydrostatic relation is
δη φ = −αd µd .
With the ARW vertical coordinate η, the model lid pt , and the column dry pressure known at
each (i, j, k) location, the 3-dimensional arrays are interpolated.
47
Vertical calculations are always interpolations in the free atmosphere, up to the model lid;
however, near the model surface, it is possible to have an inconsistency between the input
surface pressure (which is largely based on the input surface elevation) and the ARW surface
pressure (which could have a significantly higher-resolution topography). These inconsistencies
may lead to an extrapolation. The default behavior for extrapolating the horizontal winds and
the relative humidity below the known surface keeps the values constant, with zero vertical
gradient. For potential temperature, a default of 6.5 K/km lapse rate is applied. Vertical
interpolation of the geopotential field is optional and is handled separately. Since a known lower
boundary condition exists (the geopotential is defined as zero at the pressure at sea-level), no
extrapolation is required.
µ0d = µd − µd , (5.9)
Starting with the reference state fields (5.4, 5.6, and 5.7) and the perturbation (5.9), the pertur-
bation fields for pressure and inverse density are diagnosed. The pressure perturbation includes
moisture, and is diagnosed from the hydrostatic equation
0 0
δη p = µ d 1 + q v + q v η µ d ,
η
which is integrated down from the model top (where p0 = 0) to recover p0 . The total dry inverse
density is given as
0 − CCv
Rd pd + pd p
αd = θm ,
p0 p0
which defines the perturbation field for inverse density
αd0 = αd − αd .
The perturbation geopotential is diagnosed from the hydrostatic relation
δη φ0 = − µd αd0 + µ0d αd
by upward integration, using the terrain elevation as the lower boundary condition.
48
For real-data cases, the specified lateral boundary condition for the coarse grid is supplied
by an external file that is generated by program real. This file contains records for the fields u,
v, θ, qv , φ0 , and µ0d that are used by ARW to constrain the lateral boundaries (other fields are
in the boundary file, but are not used). The lateral boundary file holds one less time period
than was processed by WPS. Each of these variables has both a valid value at the initial time of
the lateral boundary time, and a tendency term to get to the next boundary time period. For
example, assuming a 3-hourly availability of data from WPS, the first time period of the lateral
boundary file for u would contain data for both coupled u (map scale factor and µd interpolated
to the variable’s staggering) at the 0 h time
µd x u
U0h = ,
mx 0h
which would take a grid point from the initial value to the value at the 3-hour simulation time.
The horizontal momentum fields are coupled both with µd and the inverse map factor. The
other 3-dimensional fields (θ, qv , and φ0 ) are coupled only with µd . The µ0d lateral boundary
field is not coupled.
Each lateral boundary field is defined along the four sides of the rectangular grid (loosely
referred to as the north, south, east, and west sides). Boundary values and tendencies for vertical
velocity and the non-vapor moisture species are included in the external lateral boundary file,
but act as place-holders for the nested boundary data for the fine grids. The width of the lateral
boundary zone along each of the four sides is user-selectable at run-time (as shown in Fig. 6.1).
49
5.3 Digital Filtering Initialization
ARW provides a digital filtering initialization (DFI) to remove noise from the model initial
state. This noise results from imbalances between the interpolated mass and wind fields, and
the differences between the coarser topography of the first-guess data set and relatively higher
resolution topography of ARW. DFI is applied to output from the real preprocessor before the
model simulation begins. If data assimilation is performed using WRF-Var, DFI is applied to
analysis produced by the WRF-Var system, rather than the output of program real.
Under the assumption that any noise is of higher frequency than meteorologically significant
modes, DFI attempts to remove this noise by filtering all oscillations above a specified cutoff
frequency. Accordingly, filters in the ARW DFI are low-pass digital filters, which are applied
to time-series of model fields. The initialized model state is the output of the filter at some
prescribed time, (e.g., the analysis time). Time series of model states are generated through
combinations of integration types that may include adiabatic, backward integration and diabatic,
and/or forward integration in the model, along with the chosen DFI scheme (which determines
the specific combination of integration). Three DFI schemes are available: digital filter launch
(DFL; Lynch and Huang (1994)), diabatic DFI (DDFI; Huang and Lynch (1993)), and twice
DFI (TDFI; Lynch and Huang (1994)).
where yn is the output of the filter at time n, the hk are the coefficients of the filter, and
{. . . , xn−1 , xn , xn+1 , . . .} is the sequence of input values to be filtered; such a filter is said to have
span 2N + 1.
One method for deriving the coefficients of a nonrecursive digital filter is the windowed-
sinc method, described in the context of DFI by Lynch and Huang (1992). In the ARW DFI,
either the Lanczos, Hamming, Blackman, Kaiser, Potter, or Dolph-Chebyshev windows may
be used; the Dolph-Chebyshev window is described by Lynch (1997). However, when a filter
with a shorter span is desired, another nonrecursive digital filter - the Dolph filter - may be
used. Lynch (1997) describes the construction of the Dolph filter, and demonstrates that this
filter has properties nearly indistinguishable from those of an optimal filter, which minimizes
the maximum difference between a filter’s transfer function and an ideal transfer function in the
pass and stop bands.
The only recursive filter in ARW DFI is the second-order Quick-Start filter of Lynch and
Huang (1994). In general, a recursive digital filter that depends only on past and present values
of input, and on past values of output, is of the form
50
Digital Filter Launch
t=0
Figure 5.2: An illustration showing the three available DFI schemes: digital filter launch, dia-
batic digital filter initialization, and twice digital filter initialization.
N
X N
X
yn = hk xn−k + bk yn−k . (5.11)
k=0 k=1
However, this form is inconvenient when inputs to the filter consist of model states, and storing
many such states should be avoided. Lynch and Huang (1994) show how this type of recursive
filter can be reformulated to the same form as a nonrecursive filter, and thus, the second-order
Quick-Start filter follows the same form as (5.10).
DFL
In the DFL scheme, forward integration with full model physics and diffusion begins at the
initial time and continues for 2N time steps, during which time a filtered model state valid N
time steps beyond the analysis time is computed as in (5.10). Then, the initialized simulation
is launched from the midpoint of the filtering period. For any model state X, let [X]D n give
the model state after diabatically integrating n time steps forward in time; we emphasize that
the superscript D indicates diabatic integration, in contrast to adiabatic integration. Then, the
DFL scheme is expressed as
2N
X
Xini = hn [Xana ]D
n , (5.12)
n=0
where Xini is the initialized model state, Xana is the model analysis or model initial state
generated by the real preprocessor, and hn represents filter coefficients.
51
DDFI
To produce an initialized state valid at the model analysis time, the DDFI scheme begins with
an adiabatic, backward integration for N time steps, followed by a diabatic, forward integration
for 2N time steps, during which filtering takes place. This filtered state is valid at the model
analysis time. Letting [Xana ]A
−n denote the model state after adiabatic, backward integration for
n time steps from the model analysis or model initial state, Xana , the DDFI scheme is expressed
as
2N
X h iD
Xini = hn [Xana ]A
−n , (5.13)
n
n=0
where Xini is the initialized model state valid at the model analysis time.
TDFI
The TDFI scheme involves two applications of the digital filter. Adiabatic, backward integration
proceeds from the model initial time for 2N time steps, during which a filter is applied. The
filtered state is valid at time −N ∆t. From this filtered state, a forward, diabatic integration is
launched. The second integration proceeds for 2N time steps, during which a second filter is
applied, giving a filtered model state valid at this model analysis time. The TDFI scheme is
expressed as
2N
" 2N #D
X X
Xini = hn hn0 [Xana ]A
−n . (5.14)
n=0 n0 =0 n
th
3rd order: (q xadv )i−1/2 = (q xadv )4i−1/2
1
− sign(U ) (qi+1 − qi−2 ) − 3(qi − qi−1 )
12
th
5th order: (q xadv )i−1/2 = (q xadv )6i−1/2
1
+ sign(U ) (qi+2 − qi−3 ) − 5(qi+1 − qi−2 ) + 10(qi − qi−1 ) .
60
When specified boundary conditions are used, as in Section 6.4, the model boundaries before
the model initial time are taken to be the same as those valid at the model initial time. We note
that, with a negated time step, the linear ramping functions F1 and F2 of (6.1) change sign,
and, consequently, so does the sign of the tendency for a prognostic variable ψ.
52
Chapter 6
ARW provides several lateral boundary conditions suitable for idealized flow, in addition to a
specified lateral boundary condition for real-data simulations. These choices are handled via a
run-time option in the Fortran namelist file. For nesting, all fine domains use the nest time-
dependent lateral boundary condition where the outer row and column of the fine grid are
specified from the parent domain, described in Section 7.3. The remaining lateral boundary
options are exclusively for use by the outer-most parent domain.
for each integer (n, m). The periodicity lengths (Lx , Ly ) are
[(dimension of the domain in x) - 1]∆x and [(dimension of the domain in y) - 1]∆y.
δτ U 00 = −U ∗ δx u,
53
horizontal difference operator δx is evaluated in a one-sided manner using the difference between
the boundary value and the value one grid-point into the grid, from the boundary. cb is the
phase speed of the gravity waves that are to be radiated; it is specified as a model constant (for
more details see Klemp and Lilly, 1978; Klemp and Wilhelmson, 1978).
For scalars and non-normal momentum variables, the boundary-perpendicular flux diver-
gence term is replaced with
δx (U ψ) = U ∗ δx ψ + ψ δx U,
where U ∗ = min(U, 0) at the x = 0 + ∆x/2 (western) scalar boundary, U ∗ = max(U, 0) at the
x = L − ∆x/2 (eastern) boundary, and likewise for the south and north boundaries, using V . As
was the case for the momentum equations, the horizontal difference operator δx is evaluated in a
one-sided manner using the difference between the boundary value and the value one grid-point
into the grid from the boundary.
where xb is the location of the symmetry boundary. All other variables satisfy the relation
54
Real-Data Lateral Boundary Condition: Location of Specified and Relaxation Zones
Specified
Columns
Relaxation
North Columns
N=1
N=2
Specified N=3
Rows N=4
N=5
West East
1 2 3 4 5 5 4 3 2 1
Relaxation
Rows
N=5
N=4
N=3
N=2
N=1
South
Figure 6.1: Specified and relaxation zones for a grid with a single-specified row and column,
plus four rows and columns for the relaxation zone. These are typical values used for a specified
lateral boundary condition for a real-data case. In this figure, the weight of the relaxation term
would be identically zero for the fifth row or column in from the boundary edge.
along the outer edge of the outer-most grid is entirely specified by temporal interpolation, using
data from an external model). The second region of the lateral boundary for the coarse grid
is the relaxation zone, which is where the model is nudged or relaxed towards the large-scale
forecast (e.g., rows and columns 2 through 5 in Fig. 6.1). The size of the relaxation zone is a
run-time option.
The specified lateral boundary condition for the coarse grid requires an external file, gener-
ated during the same pre-processing as the initial condition file. Let ψ be any prognostic value
having a lateral boundary entry, after Davies and Turner (1977),
∂t ψ n
= F1 (ψLS − ψ) − F2 ∆2 (ψLS − ψ), (6.1)
where n is the number of grid points in from the outer row or column along the boundary
(SpecZone + 1 ≤ n ≤ SpecZone + RelaxZone − 1; see Fig. 6.1) and ψLS is the large-scale value
obtained by spatial and temporal interpolation from the external analysis or model forecast by
WPS . ∆2 is a 5-point horizontal smoother applied along η-surfaces. The weighting function
coefficients F1 and F2 in (6.1) are given by
1 SpecZone + RelaxZone − n
F1 = ,
10∆t RelaxZone − 1
1 SpecZone + RelaxZone − n
F2 = ,
50∆t RelaxZone − 1
where n extends only through the relaxation zone (SpecZone+1 ≤ n ≤ SpecZone+RelaxZone−
1). F1 and F2 are linear ramping functions with a maximum at the first relaxation row or column
55
nearest to the coarse-grid boundary (just inside the specified zone). These linear functions can
optionally be multiplied by an exponential function for smoother behavior when broad boundary
zones are needed. This is achieved through an extra factor.
1 SpecZone + RelaxZone − n
F1 = exp[−SpecExp(n − SpecZone − 1)]
10∆t RelaxZone − 1
1 SpecZone + RelaxZone − n
F2 = exp[−SpecExp(n − SpecZone − 1)]
50∆t RelaxZone − 1
where SpecExp is an inverse length scale in grid lengths.
On the coarse grid, the specified boundary condition applies to the horizontal wind com-
ponents, potential temperature, φ0 , µ0d , and water vapor. The lateral boundary file contains
enough information to update boundary zone values through the entire simulation period. The
momentum fields are coupled with µd and the inverse map factor (both at the specific variable’s
horizontal staggering location), and the other 3-dimensional fields are coupled with µd . The µ0d
field is not coupled in the lateral boundary file.
Several fields within the model have specified lateral boundary conditions that are associated
particularly with a known variable. Vertical velocity has a zero gradient boundary condition
applied in the specified zone (the outer-most row and column). Users are able to control the
behavior of the treatment of the lateral boundary conditions for scalar variables that are subject
to transport. By default, microphysical variables (except vapor) and all other scalars (such
as chemistry constituents, passive tracers, scheme-specific aerosols, and number concentration
fields) have flow-dependent boundary conditions applied in the specified zone. This boundary
condition is defined as zero-value on inflow and zero-gradient on outflow. Since these boundary
conditions require only information from the interior of the grid, these scalar variables are not
required to be in the lateral boundary condition file. If sufficient scalar data is provided to the
lateral boundary condition file, the same technique of relaxation towards temporally-interpolated
variables is employed, identically as with the dynamic variables.
56
Chapter 7
Nesting
ARW provides the capability to focus the area of a simulation via nesting options. ARW
supports horizontal nesting that allows resolution to be enhanced over a region of interest by
introducing an additional grid (or grids) into the simulation. This subdomain with increased
horizontal resolution is completely contained within the parent domain (this is also commonly
referred to as the coarse domain and abbreviated CG). The subdomain is commonly referred
to as the child domain or the fine grid (FG). Vertical refinement options are available, where a
child domain may have a different vertical structure than the parent domain. Parent and child
domains always extend from the lower model surface, through to the upper model lid. In this
way, the vertical capability is more accurately referred to as a refinement, as opposed to nesting.
In practice, both of these terms will be used interchangeably when referring to a child domain
that has enhanced vertical resolution, compared to the parent domain.
As with all domains within ARW, the horizontally or vertically nested grids (child domains,
fine-grid domains) are rectangular and are aligned with the parent (coarser) grid within which
they are nested. The horizontal nested grids allow any integer spatial (∆xcoarse /∆xf ine ) and
temporal refinements of the parent grid (the spatial and temporal refinements are usually, but
not necessarily the same). The vertical refinement capability offers two options: either an inte-
ger factor of the number of levels in the parent to define the child domain’s vertical structure
(for example, a doubling or tripling of the vertical resolution), or a completely independent set
of vertical levels prescribed for the parent and the child domain. The horizontal nesting capa-
bility is, in many ways, similar to implementations provided in other mesoscale and cloud-scale
models (e.g. MM5, ARPS, COAMPS). Vertical refinement options are described in Mahalov
and Moustaoui (2009) and Daniels et al. (2016).
In this chapter we describe the various horizontal and vertical nesting options available in
ARW and the exchange of data between the grids.
57
1-Way ARW Simulation 2-Way ARW Simulation
Run coarse grid (CG) Both CG and FG simulations Both CG and FG simulations
simulation run within the same WRF run within the same WRF
execution execution
Recursively feedback FG to
CG
Figure 7.1: 1-way (parent provides boundary data to child) and 2-way (parent provides bound-
ary data to child, then child feeds information back to the parent) nesting options in ARW.
nest, the only information exchange is from the coarse grid to the fine grid; hence, the name
1-way nesting. In 2-way nest integration, the fine-grid solution replaces the coarse-grid solution
for coarse grid points that lie inside the fine grid. This information exchange between the
grids is now in both directions (coarse-to-fine for the fine-grid lateral boundary computation,
and fine-to-coarse during the feedback at each coarse-grid time step); hence, the name 2-way
nesting.
The 1-way nest set-up may be run in one of two different methods. One option is to produce
the nested simulation as two separate ARW simulations, as described in the leftmost box in
Fig. 7.1. In this mode, the coarse grid is integrated first and once the coarse-grid time step is
completed, output from the coarse-grid integration is processed to provide boundary conditions
for the nested run (usually at a much lower temporal frequency than the coarse-grid time step).
This is followed by the complete time integration of the fine (nested) grid. Hence, this 1-way
option is equivalent to running two separate simulations with a processing step in between. Also,
with separate grid simulations, an intermediate re-analysis (such as via WRFDA, see Section
11) may be included.
The second 1-way option (lockstep, with no feedback), depicted in the middle box in Fig.
7.1, is run as a traditional simulation with two (or more) grids integrating concurrently, except
with the feedback runtime option turned off. This option provides lateral boundary conditions
to the fine grid at each coarse-grid time step, which is an advantage of the concurrent 1-way
method (no feedback).
58
Fine Grid Initialization Options
ARW supports several strategies to horizontally refine a coarse-grid simulation with the intro-
duction of a nested grid. When using concurrent 1-way and 2-way nesting, several options for
initializing the fine grid are provided.
• All fine-grid variables (both meteorological and terrestrial) are interpolated from the coarse
grid, which is useful when a fine grid starts later in the coarse-grid forecast.
• All fine-grid variables are input from an external file that has high-resolution information
for both the meteorological and terrestrial fields. This is a standard set-up when the
fine-grid terrestrial fields are expected to impact the forecast.
• The fine grid may have some of the variables initialized with a high-resolution external
data set, while other variables are interpolated from the coarse grid (for example, this
would permit the improved analysis from the WRFDA initialization of the coarse grid’s
meteorological fields to remain consistent with the fine grid). This option allows the fine
grid to start later in the coarse grid’s forecast, but with the advantage of higher-resolution
static fields.
• For a moving nest, external orography and landuse files may be used to update the fine-grid
domain as it moves over land.
These fine-grid initialization settings are user-specified at run-time, and ARW allows nested grids
to instantiate and cease during any time that the fine grid’s parent is still integrating. While
cost savings are evident when starting the fine-grid domain at a later time, that advantage must
be weighed against the impact of relatively coarse and inconsistent fields for both masked and
meteorological variables on the fine grid.
59
1 (a) (b)
1
2 3 4 2 3
(c) (d)
1
2 x 1
x
2 4 3
3
Figure 7.2: Various nest configurations for multiple grids. (a) Telescoping nests. (b) Nests at
the same level with respect to a parent grid. (c) Overlapping grids: not allowed with feedback
activated. (d) Inner-most grid has more than one parent grid: not allowed
refinement ratio), though the model does allow the time step refinement ratio to differ from the
spatial refinement ratio. Also, nested grids on the same level (i.e., children who have the same
parent) may have different spatial and temporal refinement ratios. For example, in Fig. 7.2b,
the horizontal grid resolution for domain 1 may be 90 km, while domain 2 may be 45 km, and
domain 3 may be 30 km.
Moving Nests
The moving nest capabilities in ARW are simply extensions to the suite of nesting options. All
descriptions covering the specifics for fine-grid domains (initialization, feedback, configurations,
staggering, lateral boundaries, etc.) also apply to moving nests. In general, all nests in an ARW
forecast are eligible to be moving nests. ARW provides two methods to allow nests to move
during model integration: specified and automatic. For both types of moving nests, multiple
levels of domains may move.
For a specified move, the timing of a nest move and the extent of each lateral move is defined
entirely by the user through the namelist. This manual process is tedious and cumbersome, and
is rarely used in practice. This capability was required to eventually develop a fully functioning
moving nest capability that automatically tracks a cyclone.
For the automatically moving nest, the fine grid is initialized to cover a well-defined vortex,
and the nest moves to maintain this vortex in the center of the fine grid. The fine grid follows a
height minimum, and is constrained so as to not too closely approach the parent boundary. This
option provides a substantial cost savings when the area of horizontal refinement is reduced in
60
size to cover only the physical extent of a cyclone, and not the much larger domain necessary
to contain a moving system. The typical instance for utilizing a moving nest is during tropical
cyclone tracking via the automatic vortex-following technique. With this automated cyclone
tracking, the inner-most domain is responsible for steering the movement of the coarser, moving
grids.
After a nested domain has moved a parent grid-cell distance, the majority of the fine-grid
data in the domain is still valid. Data that are not along the outer row or column of the nested
domain is shifted to the new location in that domain; not interpolated from the parent domain.
Once a domain moves, data in the outer row or column fall into one of two categories: discarded
data on the trailing edge, or horizontally interpolated data on the leading edge in the direction
of the nest move. Only elevation of the topography and landuse categories are eligible to be
used as input from high-resolution files, instead of being interpolated from the parent domain.
61
V V V V
U q U q U q U q U
V V V V V V V V
U q U q U q U q U q U q U
V V V V V V
U q U q U q U q U q U q U q U q U
V V V V V V
U q U q U q U q U q U q U
V V V V V V V V
U q U q U q U q U q U q U
V V V V V V
U q U q U q U q U q U q U q U q U
V V V V V V
U q U q U q U q U q U q U
V V V V V V V V
U q U q U q U q U
V V V V
Figure 7.3: Arakawa-C grid staggering for a portion of a parent domain, with an imbedded nest
domain using a 3:1 grid size ratio. Solid lines denote coarse grid cell boundaries, and dashed
lines are the boundaries for each fine grid cell. The horizontal components of velocity (“U” and
“V”) are defined along the normal cell face, and the thermodynamic variables (“θ”) are defined
at the center of the grid cell (each square). The bold typeface variables along the interface,
between the coarse and the fine grid, define the locations where the specified lateral boundaries
for the nest are in effect.
62
V V V V
U q U q U q U q U
V V V V V V
U q U q U q U q U
U q V V V V q U
U q U q U q U q U
V V V V V V
U q U q U q U q U
U q V V V V q U
U q U q U q U q U
V V V V V V
U q U q U q U q U
V V V V
Figure 7.4: Similar to Fig. 7.3, but with a 2:1 grid-distance ratio.
the coarse grid. These fields include most 3-dimensional and 2-dimensional arrays. For the
horizontal momentum fields averaged back to the coarse grid in the feedback, the mean of three
(for example, due to the 3:1 grid-distance ratio in the example shown in Fig. 7.3) fine grid points
is fed back to the coarse grid from along the coincident cell face. When using odd ratios, fields
that are masked due to the land/sea category are fed back directly from the coincident points.
Masked fields include soil temperature and sea ice. It is not reasonable to average neighboring
locations of soil temperature on the fine grid if the coarse grid point to which values are being fed
back is a water value. Similarly, averaging several sea ice values on the fine grid does not make
sense if some of the neighboring points included in the mean are fine grid land points. Only
masked fields use the feedback method in which a single point (such as for land use category) or
a mean of valid points (such as for soil temperature) from the fine grid is assigned to the coarse
grid.
One difference between odd and even grid-distance ratios is in the feedback from the fine
grid to the coarse grid. No coincident points exist for the single point feedback mechanisms for
even grid distance ratios (such as is used for the land/sea masked 2D fields). For a 2:1 even
grid distance ratio, Figure 7.4 shows that each coarse grid point has four fine grid cells that
are equally close, and therefore four equally-eligible grid points for use as the single fine-grid
point that feeds back to the coarse grid. The single-point feedback is arbitrarily chosen as the
south-west corner of the four neighboring points. This arbitrary assignment to masked fields
implies that even grid distance ratios are better suited for idealized simulations where masked
fields are less important.
63
7.3 Nested Lateral Boundary Conditions
For the fine grid with 2-way nesting or 1-way nesting (using a concurrent ARW simulation, see
Fig. 7.1), boundary conditions are specified by the parent grid at every coarse-grid time step.
The nest lateral boundary condition behaves similarly to the specified boundary condition for
real-data cases (see Section 6.4), but the relaxation zone is not active. Prognostic variables
are entirely specified in the outer row and column of the fine grid through spatial and tem-
poral interpolation from the coarse grid (the coarse grid is stepped forward in time, prior to
advancement of any child grid of that parent).
Nest Instantiation
The fine grid is instantiated as a child of a parent grid at the requested start time. This
initialization is within the integration step for the parent grid, meaning that no child grid
integration can begin if the parent is not active. To fill in the correct meteorological fields, a
default initialization routine is called to horizontally interpolate coarse-grid data to the fine-grid
locations using a monotone interpolation scheme (described in Smolarkiewicz and Grell, 1992)
for most fields (i.e., the same scheme employed for generating the fine grid lateral boundary
conditions) and a simple linear interpolation, or averaging scheme, for masked or categorical
fields. For fields that are masked with the land/sea background, such as land-only fields (e.g.,
snow), or water-only fields (e.g., sea ice), the interpolator needs to know what field defines the
template for masking (such as for the land use category). Part of the automatic code generation
handles calling each field with its associated interpolator.
64
Integrate parent grid one time step
If nest grid start time
(1) Horizontally interpolate parent to child grid
(2) Optionally input high-resolution child data
(3) Compute child reference state
(4) Feedback Child Initial Data to Parent Grid
(5) Recompute parent reference state
End If nest grid start time
Solve time step for parent grid (see Fig. 3.1)
If nest grid move time and FG away from CG boundary
(1) Move nest grid (vortex following or prescribed)
(2) Horizontally translate data due to grid shift
(3) Horizontally interpolate parent to child grid (along new boundary)
(4) Optionally input high-resolution child topo-landuse data
(5) Compute child reference state
(6) Feedback child initial data to parent grid
(7) Recompute parent reference state
End If nest grid move time
While existing nest grids to integrate
(1) Lateral forcing from parent grid to child
(2) Integrate child grid to current time of parent grid
(3) Feedback child grid information to parent grid
End While existing nest grids to integrate
End Grid Integrate
65
Figure 7.6: Zones of topographic blending for a fine grid. In the fine grid, the first zone is
entirely interpolated from the coarse-grid topography. In the second zone, the topography is
linearly weighted between the coarse grid and the fine grid.
where row=1 is the first row in the second zone, and the row or column nearest the outer edge
takes precedence in ambiguous corner zones. The reference variables computed from the topog-
raphy, µd and φ, are similarly-treated. Blended arrays are required to compute the reference
state for the fine grid. The blending along the inner rows and columns ramps the coarse-grid
reference state to the fine-grid reference state for a smooth transition between grids.
Feedback
So that the coarse grid and fine grid are consistent at coincident points, the fine grid values are
fed back to the coarse grid. There are two available options for feedback: either the mean of all
fine-grid cells contained within each coarse-grid cell (or cell faces in the case of the horizontal
momentum fields) is fed back, or a single-point feedback is selected for masked or categorical
fields.
Subsequent to the feedback step, the coarse grid may be optionally smoothed in the area
of the fine grid. Two smoothers are available: a 5-point 1-2-1 smoother, and a smoother-
desmoother with a similar stencil size. Both the feedback and the smoothers are run one row
and column in from the interface row and column of the coarse grid (the coarse grid provides
the lateral boundary conditions to the fine grid, so the outer-most row cannot be modified).
66
Reference State
When the nest is instantiated, the initial feedback ensures that the coarse grid is consistent with
the fine grid, particularly with regards to topography and the reference state fields inside the
blended region, and for such terrestrial features as coasts, lakes, and islands. The adjustment
of elevation in the coarse grid forces a base state recalculation. The fine grid needs an initial
base state calculation, and after the terrain feedback, the coarse grid is also in need of a base
state recalculation. Note that with horizontal interpolation of the coarse grid to the fine grid,
and feedback of the fine grid to the coarse grid, the coarse-grid base state is recomputed even
without a separate fine-grid initial data file, since the coarse-grid topography is adjusted.
With completed base state computations, which follow similarly to that described for the
real-data initialization in section 5.2.2, routines return back to the integration step for the
coarse and fine grids. The fine-grid data are now properly initialized for integration and can be
advanced forward one time step.
Integration
The integration by grid is recursive. At the end of each grid’s time step, a check is made to
determine if a child grid exists for that parent and if the current time is bracketed by the child’s
start/end time. This is shown in Fig. 7.5. The integration process for the nest (step 2 under the
while loop) recursively calls the top step in the overall sequence as a parent grid itself. This is a
“depth first” traversal of the tree of grids. If a child grid does exist, that child grid is integrated
up through the current time of the parent grid.
Interpolation Options
Fields from the parent domain (coarse grid, CG) are interpolated to the child domain (fine
grid, FG) prior to the CG data being used by the FG. There are several occurrences of this
interpolation.
• At initialization, or upon a domain move, horizontally interpolate all non-masked fields
from CG to FG.
• At initialization or upon a domain move, topography (and fields derived from topography)
within the FG domain are blended with data from the CG, as shown in Fig. 7.6.
• At the completion of each CG timestep, CG fields are interpolated for use by the FG.
This is primarily for the construction of lateral boundary conditions for the FG domain;
however, optional interpolations are available for the fields that are only computed on the
coarsest grid (examples include stochastic forcing and simple ocean models that are only
defined on the outer-most domain).
ARW supports a number of horizontal interpolation options, available via the user-defined
namelist.
67
68
Chapter 8
Physics
This chapter outlines the physics options available in the ARW. The WRF physics options fall
into several categories, each containing several choices. The physics categories are (1) micro-
physics, (2) cumulus parameterization, (3) planetary boundary layer (PBL), (4) land-surface
model, and (5) radiation. Diffusion, which may also be considered part of the physics, is de-
scribed in Chapter 4.
The physics section is separated from the dynamics in the solver by the use of physics drivers.
These drivers are between dynamical-core-dependent parts of the solver: a pre-physics prepa-
ration of variables and post-physics modifications of the tendencies. The physics preparation
involves filling arrays with physics-required variables that include the temperature, pressure,
heights, layer thicknesses, and other state variables in MKS units at half-level grid points and
on full levels. The velocities are also de-staggered so that the physics part is independent of
the dynamical solver’s velocity staggering. Physics packages compute tendencies for the velocity
components (un-staggered), potential temperature, and moisture fields. A post-physics step will
re-stagger these tendencies as necessary, couple tendencies with coordinate metrics, and convert
to variables or units appropriate to the dynamics solver.
In the first Runge-Kutta step, prior to the acoustic steps (see Fig. 3.1, step(1)), tendencies
are computed for radiation, surface, PBL, and cumulus physics. These tendencies are then held
fixed through the Runge-Kutta steps. Microphysics is computed after the last Runge-Kutta
step (see Fig. 3.1, step (9)) in order to ensure proper saturation conditions at the end of the
time-step.
The initialization of the physics is called prior to the first model step. This initialization
may include reading in data files for physics tables or calculating look-up tables of functions.
Each physics module includes an initialization routine for this purpose. Often physics packages
will have many of their own constants that are included in their own module, while common
physical constants are passed in from the physics drivers.
8.1 Microphysics
Microphysics includes explicitly resolved water vapor, cloud, and precipitation processes. The
model is general enough to accommodate any number of mass mixing-ratio variables, and other
quantities such as number of particles per unit dry air mass. Four-dimensional arrays with three
spatial indices and one species index are used to carry such scalars. Memory, i.e., the size of the
69
fourth dimension in these arrays, is allocated depending on the needs of the scheme chosen, and
advection of the species also applies to all those required by the microphysics option. In the
current version of ARW, microphysics is carried out at the end of the time-step as an adjustment
process, and so does not provide tendencies but directly updates the state variables. The ratio-
nale for this is that condensation adjustment should be at the end of the time-step to guarantee
that the final saturation balance is accurate for the updated temperature and moisture. How-
ever, it is also important to have the latent heating forcing for potential temperature during the
dynamical sub-steps, and this is done by saving the microphysical heating as an approximation
for the next time-step as described in Section 3.1.4.
Currently, the sedimentation process is accounted for inside the individual microphysics
modules, and, to prevent instability in the calculation of the vertical flux of precipitation, either
a smaller time step is allowed or a Lagrangian scheme is used. The saturation adjustment is
also included inside the microphysics. In the future, however, it might be separated into an
individual subroutine to enable the remaining microphysics to be called less frequently than the
model’s advection step for efficiency.
The various microphysics options have differing numbers of moisture variables, depending on
the ice-phase and mixed-phase processes included. Mixed-phase processes are those that result
from the interaction of ice and water particles, such as riming that produces graupel or hail.
As a general rule, for grid sizes less than 10 km, where updrafts may be resolved, mixed-phase
schemes should be used, particularly in convective or icing situations. For coarser grids it may
not be worth the added expense of these schemes because riming is not likely to occur with
the relatively weak vertical motion that is resolved. Many schemes are also double-moment for
some species and include number per unit dry air mass for those as extra advected variables.
70
proaches is that a diagnostic relation is used for ice number that is based on ice mass content
rather than temperature. The computational procedures are described in Hong and Lim (2006).
As with WSM5 and WSM6, the fall terms are computed using a Lagrangian technique instead
of time-splitting that was used in earlier versions of these schemes. The WSM3 scheme predicts
three categories of moist variables: water vapor, cloud water/ice, and rain/snow, making it a
so-called simple-ice scheme. It follows Dudhia (1989) in assuming cloud water and rain for tem-
peratures above freezing, and cloud ice and snow for temperatures below freezing. This scheme
is computationally efficient for the inclusion of ice processes, but lacks supercooled water and
gradual melting rates.
71
8.1.7 WSM7 and WDM7 schemes
The extension of WSM6 and WDM6 to include hail as a separate category from graupel is
also available as a choice in WRF as WSM7 and WDM7 respectively (Bae et al., 2018). These
schemes, by allowing for hail, also can give more intense precipitation rates. As with the other
WDM schemes, WDM7 only affects the cloud and rain water processes.
72
First, there is an option to choose either graupel or hail as the third class of ice in the so-
called 3ICE scheme (McCumber et al., 1991). Graupel has a relatively low density and a high
intercept value (i.e., more numerous small particles). In contrast, hail has a relative high density
and a low intercept value (i.e., more numerous large particles). These differences can affect not
only the description of the hydrometeor population and formation of the anvil-stratiform region
but also the relative importance of the microphysical-dynamical-radiative processes.
Second, new saturation techniques (Tao et al., 1989, 2003) were added. These saturation
techniques are basically designed to ensure that super saturation (sub-saturation) cannot exist
at a grid point that is clear (cloudy).
Third, all microphysical processes that do not involve melting, evaporation or sublimation
(i.e., transfer rates from one type of hydrometeor to another) are calculated based on one
thermodynamic state. This ensures that all of these processes are treated equally.
Fourth, the sum of all sink processes associated with one species will not exceed its mass.
This ensures that the water budget will be balanced in the microphysical calculations .
The Goddard microphysics has a simpler option, which is equivalent to a two-ice (2ICE)
scheme having only cloud ice and snow. This option may be sufficient for coarse resolution
simulations (i.e., > 10 km grid size). The two-class ice scheme could be applied for winter and
frontal convection.
73
8.1.13 Milbrandt-Yau Double-Moment scheme
The Milbrandt and Yau (2005) scheme is a fully double-moment scheme that carries graupel
and hail as separate species. It therefore has 12 prognostic variables in addition to water vapor.
74
8.1.19 Jensen ISHMAEL microphysics
The ISHMAEL (Ice-Spheroids Habit Model with Aspect-ratio EvoLution) scheme (Jensen et al.,
2017) can predict ice-particle habits that are defined by their aspect ratios and volumes and how
these develop through deposition and riming. Up to three independent habits are carried which
contain generalized particles that represent the whole evolution from ice to snow to graupel
together with their numbers per unit dry air mass.
75
• Downdraft changes:
– Source layer is the entire 150 – 200 mb deep layer just above cloud base.
– Mass flux is specified as a fraction of updraft mass flux at cloud base. Fraction is
a function of source layer RH rather than wind shear or other parameters, i.e., old
precipitation efficiency relationship not used.
– Detrainment is specified to occur in updraft source layer and below.
• A new way to compute perturbation temperature added to the testing updraft parcel is
implemented based on Ma and Tan (2009). This perturbation temperature is a function
of horizontal and vertical moisture advection, rather than vertical velocity. This changes
the trigger function in the original scheme to work better in a weakly forced environment.
• Subgrid cloud fraction due to both deep and shallow convection is estimated and passed
to the radiation schemes (Alapaty et al., 2012).
76
many details and/or parameter values differ from those recommended in Betts (1986) and Betts
and Miller (1986). Recently, attempts have been made to refine the scheme for higher horizontal
resolutions, primarily through modifications of the triggering mechanism. In particular:
• A floor value for the entropy change in the cloud is set up below which the deep convection
is not triggered;
• In searching for the cloud top, the ascending parcel mixes with the environment; and
• The work of the buoyancy force on the ascending parcel is required to exceed a prescribed
positive threshold.
Grell-3 Scheme
The Grell-3 scheme was first introduced in Version 3.0. It shares a lot in common with the
Grell-Devenyi in scheme, being based on an ensemble mean approach, but the quasi-equilibrium
approach is no longer included among the ensemble members. The scheme is distinguished from
other cumulus schemes by allowing subsidence effects to be spread to neighboring grid columns,
making the method more suitable to grid sizes less than 10 km, while it can also be used at
larger grid sizes where subsidence occurs within the same grid column as the updraft.
Grell-Freitas Scheme
The Grell-Freitas scheme (Grell and Freitas, 2014) is based on Grell-Devenyi scheme which
considers a stochastic approach to cumulus convection, and modified to work across grid-sizes
from mesoscale to convective scales. The scheme relates the convective updraft fraction to the
entrainment rate, which in turn is related the cloud radius. As the grid size decreases, the
fractional updraft area increases, which is equivalent to decreasing the unit mass flux required
to stabilize the atmosphere.
77
8.2.4 Simplified Arakawa-Schubert Schemes
This group of schemes is based on simplified Arakawa-Schubert scheme (SAS) (Grell, 1993) which
considers only one type of cloud instead of an ensemble. They all use the quasi-equilibrium
closure. They vary by how they handle convective downdrafts, shallow convection, momen-
tum transport (all three not included in the original AS and last two not in SAS), entrain-
ment/detrainment and scale-aware aspects.
78
8.2.5 Zhang-McFarlane Scheme
The Zhang-McFarlane scheme is a mass-flux scheme and it is only for deep convection (Zhang and
McFarlane, 1995). It follows the ideas of Arakawa-Schubert, but makes simplified assumptions.
One of the assumptions is that instead of using a spectrum of cloud plumes, it specifies the
distribution of the updraft assuming that they all have the same cloud base mass flux, but each
has a characteristic entrainment rate. It also assumes that the convection removes CAPE at
an exponential rate with a specified adjustment time scale. The scheme detrains cloud liquid
and ice, and considers momentum transport. The scheme is ported from CAM4, and tested
with other CAM physics. For example, it can be used together with Park-Bretherton shallow
convection option and CAM microphysics.
Tiedtke Scheme
The Tiedtke scheme is a mass-flux scheme, and it parameterizes deep, shallow and midlevel
convection. It represents the cloud ensemble by a bulk cloud model, and considers entrainment
and detrainment and downdrafts. It uses a CAPE closure to determine the strength of the
deep and midlevel convection, and surface evaporation for shallow convection. In the version
implemented in WRF (Zhang and Wang, 2011), turbulent entrainment and detrainment are
added and turbulent entrainment for shallow convection is increased to promote boundary layer
cloud formation. Another modification is to change how detrained cloud at the top is treated:
in the original scheme, it evaporates immediately; and in the current scheme, it is added to the
grid scale. To limit the deep convection in drier regions, a minimum relative humidity of 80%
is imposed for the mean RH between cloud base and top.
This is an updated version of Tiedtke scheme, and it is closer to the one used in the recent
ECMWF IFS (Zhang and Wang, 2017). The updates include trigger functions for deep and
shallow convection, closures for deep and shallow convection, convective adjustment time scale,
entrainment and detrainment rates for all types of convection, conversion from cloud water/ice to
rain/snow and options for momentum transport. The entrainment for deep convection is made
to depend on environmental moisture which helps the simulated tropical systems. The updated
CAPE closure for deep convection relaxes to CAPE generated in the planetary boundary layer
which improves the diurnal precipitation over land. The convective adjustment time is no longer
a constant, but it depends on the vertical velocity averaged in the updraft and cloud depth.
79
8.3 Shallow Cumulus Parameterization
Similar to the deep cumulus parameterization, these schemes represent the sub-grid-scale trans-
port of heat and moisture in shallow, and sometimes non-precipitating, clouds. Since the scale
of shallow convective cloud is generally smaller than its deep counterpart, it should be used
in models when a deep cumulus scheme is turned off. These stand-alone shallow schemes can
therefore serve that purpose.
80
and PBL schemes. Currently, each surface layer option is tied to particular boundary-layer
options, but in the future more interchangeability and options may become available. Note that
some boundary layer schemes (ACM2 and MRF) require the thickness of the surface layer in the
model to be representative of the actual surface layer (e.g. 50-100 meters) while most others can
have thin surface layers. The schemes use Monin-Obukhov similarity theory with variations in
the stability functions and methods for computing roughness lengths. Similarity theory relates
the information at the lowest model level to the surface via a stability dependent, approximately
log, profile of wind and scalars. Diagnostic outputs of 2m and 10m quantities are consistently
computed with the profiles of these schemes.
81
8.4.4 MYNN surface layer
The MYNN surface-layer scheme includes several forms for stability functions with Dyer and
Hicks (1970) used by default. There are also several options for handling thermal roughness
lengths and the fluxes over water. It is to be used with the MYNN PBL options that is part of
the RAP/HRRR physics suite.
82
vegetation, root, and canopy effects and surface snow-cover prediction. The land-surface model
provides no tendencies, but does update the land’s state variables which include the ground
(skin) temperature, soil temperature profile, soil moisture profile, snow cover, and possibly
canopy properties. There is no horizontal interaction between neighboring points in the LSM,
so it can be regarded as a one-dimensional column model for each WRF land grid-point.
83
into account phase changes of soil water (Smirnova et al., 1997, 2000). The RUC LSM also
has a multi-layer snow model with changing snow density, refreezing liquid water percolating
through the snow pack, snow depth and temperature dependent albedo, melting algorithms
applied at both snow-atmosphere interface and snow-soil interface, and simple parameterization
of fractional snow cover with possibility of grid averaged skin temperature going above freezing.
It also includes vegetation effects and canopy water. The RUC LSM has a layer approach to the
solution of energy and moisture budgets. The layer spans the ground surface and includes half
of the first atmospheric layer and half of the top soil layer with the corresponding properties
(density, heat capacity, etc.) The residual of the incoming fluxes (net radiation, latent and
sensible heat fluxes, soil heat flux, precipitation contribution into heat storage, etc.) modifies
the heat storage of this layer. An implicit technique is applied to the solution of these equations.
Prognostic variables include soil temperature, volumetric liquid, frozen and total soil moisture
contents, surface and sub-surface runoff, canopy moisture, evapotranspiration, latent, sensible
and soil heat fluxes, heat of snow-water phase change, skin temperature, snow depth and density,
and snow temperature. This option also has a sub-grid mosaic sub-option to allow for fractional
areas of different land-use categories.
84
8.5.7 Simplified Simple Biosphere Model (SSiB)
SSiB (Xue et al., 1991; Sun and Xue, 2001), the LSM from the UCLA GCM as another climate
option for land-surface physics. It has a 2-layer bulk soil temperature model, vegetation effects
at the surface, 3 layers of soil moisture with a root zone, and a 4-layer snow treatment. It also
predicts canopy temperature and moisture.
The urban canopy model estimates the surface temperature and heat fluxes from the roof,
wall and road surface. It also calculates the momentum exchange between the urban surface
and the atmosphere. If they are available, the UCM can take three different densities of urban
development using special land-use categories. Since Version 3, an anthropogenic heating diurnal
cycle has been added as an option.
85
8.5.11 Ocean Mixed-Layer Model
This can be selected as an independent option for water surfaces, and is designed for hurricane
modeling in order to simulate the cooling of the ocean underneath hurricanes. The ocean mixed-
layer model is based on that of Pollard et al. (1973). Each column is independently coupled
to the local atmospheric column, so the model is one-dimensional. The ocean part consists of
a time-varying layer, representing the variable-depth mixed layer over a fixed layer acting as a
reservoir of cooler water with a specified thermal lapse rate. In the mixed layer, the prognostic
variables are its depth, vector horizontal current, and mean temperature taken to be the sea-
surface temperature (SST). The hurricane winds drive the current, which in turn leads to mixing
at the base of the mixed layer when the Richardson number becomes low enough. This mixing
deepens and cools the mixed layer, and hence the cooler sea-surface temperature impacts the
heat and moisture fluxes at the surface, and has a negative feedback on hurricane intensity. The
model includes Coriolis effects on the current, which are important in determining the location
of maximum cooling on the right side of the hurricane track. It also includes a mixed-layer heat
budget, but the surface fluxes and radiation have much less impact than the hurricane-induced
deep mixing on the thermal balance at the time scales considered during a forecast. The ocean
mixed-layer model is initialized using the observed SST for the mixed layer, and with a single
depth representative of known conditions in the hurricane’s vicinity that may be replaced with a
map of the mixed-layer depth, if available. The initial current is set to zero, which is a reasonable
assumption given that the hurricane-induced current is larger than pre-existing ones.
86
the fluxes are combined with those of open water. Usually the depth is considered fixed and the
fraction may be updated with the sea-surface temperature periodically during the simulation as
the model has no prognostic equation for sea-ice amounts. The sea ice in Noah and NoahMP
considers 4 layers each 1 meter thick for the sea-ice energy budget.
87
dependent on the buoyancy profile, in which the PBL top is defined at the maximum entrainment
layer (compared to the layer at which the diffusivity becomes zero). A smaller magnitude of the
counter-gradient mixing in the YSU PBL produces a well-mixed boundary-layer profile, whereas
there is a pronounced over-stable structure in the upper part of the mixed layer in the case of the
MRF PBL. Details are available in Hong et al. (2006), including the analysis of the interaction
between the boundary layer and precipitation physics.
Topographic drag effects were added as an option to this PBL scheme by Jimenez and Dudhia
(2012) and improved by Lorente-Plazas et al. (2016) which modifies the model drag according
to sub-grid variance in terrain elevation and also resolved local variability.
Top-down mixing was added as an option (Wilson and Fovell, 2018) to allow for radiative-
driven downward mixing that helps the life-cycle of stratocumulus clouds and fog.
88
including shallow convection and EDMF options (Olson et al., 2019) as well as updated options
for computing mixing length scales. The schemes can be used with the MYNN or MM5 surface-
layer schemes.
An additional option related to the MYNN PBL scheme is a wind-farm parameterization
(Fitch et al., 2012) that accounts for the additional drag and turbulence generation by wind-
farm rotors.The scheme is customizable to different rotor characteristics as a function of wind
speed.
89
8.6.9 Shin-Hong PBL
Shin and Hong (2015) have developed a scale-aware PBL option based on the YSU PBL scheme.
At larger grid sizes it resembles YSU, but as the grid size becomes much less than the PBL
depth, the nonlocal term is reduced in strength to allow the resolved scales to do a fraction of
the transport consistent with resolution.
90
diagnostic outputs of the diffuse and direct components of solar radiation at the surface. Di-
agnostics for the RRTMG and CAM3 options also include TOA/surface, longwave/shortwave,
clearsky/allsky, upward/downward accumulated fluxes for radiation budgets.
Within the atmosphere the radiation responds to model-predicted cloud and water vapor
distributions, (diagnostic) cloud fraction, as well as specified carbon dioxide, ozone, and (op-
tionally) trace gas concentrations. Some schemes can handle aerosols. All the radiation schemes
in WRF currently are column (one-dimensional) schemes, so each column is treated indepen-
dently, and the fluxes correspond to those in infinite horizontally uniform planes with cloud
fractions at leach layer, which is a good approximation if the vertical thickness of the model
layers is much less than the horizontal grid length. This assumption would become less accurate
at high horizontal resolution.
In addition to radiative transfer schemes listed below, WRF also has idealized temperature
relaxation methods for the Held-Suarez global test case and the idealized tropical cyclone case.
91
properties computed by WRF-Chem or can be specified/input in other ways. The longwave
scheme has been modified at the top-of-atmosphere to account for the significant downward flux
originating from above the model top (Cavallo et al., 2011).
92
for scattered and reflected components. Ozone is considered with several climatological one-
dimensional profiles available. The scheme can also take aerosol optical properties provided by
WRF-Chem.
93
Figure 8.1: Diagram showing interactions between various physics components.
Table 8.1: Physics Interactions. Columns correspond to model physical processes: radiation
(Rad), microphysics (MP), cumulus parameterization (CP), planetary boundary layer/vertical
diffusion (PBL), and surface physics (Sfc). Rows corresponds to model variables where i and o
indicate whether a variable is input or output (updated) by a physical process.
Atmospheric Momentum io io
State or Pot. Temp. io io io io
Tendencies Water Vapor i io io io
Cloud i io o io
Precip i io o
Surface Longwave Up i o
Fluxes Longwave Down o i
Shortwave Up i o
Shortwave Down o i
Sfc Convective Rain o i
Sfc Resolved Rain o i
Heat Flux i o
Moisture Flux i o
Surface Stress i o
94
The radiation, boundary-layer and cumulus parameterization schemes all output tendencies,
but the tendencies are not added until later in the solver where their sum contributes to the
forcing for the dynamics, so from this perspective the order of call is not important. Moreover,
these physics schemes do not have to be called at the same frequency as each other or the model
time step. When lower frequencies are used, their tendencies are kept constant between calls.
This is typically done for the radiation schemes, which are too expensive to call every timestep,
and often for the cumulus schemes, for which it may only be necessary to update the convective
columns every five minutes or so. However, the surface/boundary-layer schemes are normally
called every step in ARW because this is likely to give the best results, but the option is available
to use a less frequent call that may be useful for the more expensive CLM4 LSM.
The radiation is called first because of the required radiative fluxes that are input to the
land-surface scheme. The land-surface also requires rainfall from the microphysics and cumulus
schemes, but that is from the previous time-step since it is called before the cumulus scheme
and the microphysics is at the end of the time-step. The boundary-layer scheme is necessarily
after the land-surface scheme because it requires the heat and moisture fluxes. Some cumulus
schemes (Grell options) also use the boundary-layer tendencies.
95
96
Chapter 9
97
Figure 9.1: Perturbation patterns for three different spatial scales: a) convection-scale, b) meso-
scale, c) synoptic scale.
∂X X
= D + (1 + r) Pi . (9.1)
∂t i
Here, ∂X/∂t denotes the total tendency in variable X, D is the tendency from the dynami-
cal core, Pi the tendency from the i-th physics scheme, and r is a two-dimensional, Gaussian
distributed zero-mean random perturbation field with spatial and temporal correlations (as in
Fig.??). Depending on the implementation, the SPPT scheme can use up to three patterns with
different spatial and time scales, as well as a vertical tapering function which tapers perturba-
tions near the surface and for the highest model levels. The WRF implementation by default
uses only a single perturbation field (Berner et al., 2015), and it perturbs the tendencies from
the PBL and convection schemes, but not those from the radiation or microphysics schemes. For
multiple domains, the perturbation pattern of the parent domain is interpolated to the nested
domains.
98
drag. In the WRF implementation, no dissipation weighting is used and the perturbations are
extended to the potential temperature field (Berner et al., 2011). Wind component perturba-
tions are proportional to the square root of the kinetic-energy backscatter rate, and temperature
perturbations are proportional to the potential energy backscatter rate.
A comparison of SKEBS and SPPT shows that SKEBS introduces the most model diversity
in the free atmosphere and for dynamical variables, while SPPT is most active in regions with
large tendencies (e.g., in areas with convection and near the surface) (Berner et al., 2015).
99
9.5 Stochastic Perturbations to the Lateral Boundary
Conditions
The stochastic tendencies in WRF are typically treated as physics tendencies, and they change
the perturbed fields at each time step within the WRF domain. However, since the stochastic
perturbation field is also generated on the boundary, it can be used to perturb the lateral
boundary specified and relaxation zones.
100
Chapter 10
Nudging
101
a cold-start using just the 0 hour analysis because it gives the model a chance to spin up. In
particular, the model will have six hours to adjust to topography, and produce cloud fields by
hour 0, whereas with a cold start there would be a spin-up phase where waves are produced and
clouds are developed. This method could also be combined with 3D-Var techniques that may
provide the hour 0 analysis.
Grid nudging has been added to ARW using the same input analyses as the WPS pre-
processing systems can provide. Since it works on multiple domains in a nesting configuration,
it requires multiple time-periods of each nudged domain as input analyses. Given these analyses,
the real program produces another input file which is read by the model as nudging is performed.
This file contains the gridded analysis 3d fields of the times bracketing the current model time
as the forecast proceeds. The four nudged fields are the two horizontal wind components (u and
v), temperature, and specific humidity.
The method is implemented through an extra tendency term in the nudged variable’s equa-
tions, e.g.
∂θ
= F (θ) + Gθ Wθ (θ̂0 − θ)
∂t
where F (θ) represents the normal tendency terms due to physics, advection, etc., Gθ is a time-
scale controlling the nudging strength, and Wθ is an additional weight in time or space to limit
the nudging as described more below, while θ̂0 is the time- and space-interpolated analysis field
value towards which the nudging relaxes the solution.
Several options are available to control the nudging.
a) Nudging end-time and ramping. Nudging can be turned off during the simulation, as
in dynamical initialization. Since turning nudging off suddenly can lead to noise, there is a
capability for ramping the nudging down over a period, typically 1-2 hours to reduce the shock.
b) Strength of nudging. The timescale for nudging can be controlled individually for winds,
temperature and moisture. Typically the namelist value of 0.0003 s-1 is used, corresponding
to a timescale of about 1 hour, but this may be reduced for moisture where there may be less
confidence in the analysis versus the details in the model.
c) Nudging in the boundary layer. Sometimes, since the analysis does not resolve the diurnal
cycle, it is better not to nudge in the boundary layer to let the model PBL evolve properly,
particularly the temperature and moisture fields. Each variable can therefore be selectively not
nudged in the model boundary layer, the depth of which is given by the PBL physics.
d) Nudging at low levels. Alternatively the nudging can be deactivated for any of the variables
below a certain layer throughout the simulation. For example, the lowest ten layers can be free
of the nudging term.
e) Nudging and nesting. Each of these controls is independently set for each domain when
nesting, except for the ramping function, which has one switch for all domains.
102
10.3 Flux-Adjusting Surface Analysis Nudging
This is a surface analysis nudging method that not only nudges temperature and water mixing
ratio at the first model level, but it also corrects soil temperature and moisture based on the
differences in temperature and water vapor mixing ratio at the first model level (also refered to as
Flux-Adjusting Surface Data Assimilation System or FASDAS). The differences of temperature
and water vapor mixing ratio at the first model level are converted to surface flux adjustments
that are applied to the prognostic equations of soil temperature and soil moisture. Through
this indirect nudging, the surface temperature and water vapor mixing ratio from the model will
become closer to the analysis. This method was developed for MM5 by Alapaty et al. (2008).
N
Wa2 (i, x, y, z, t)) [ao (i) − am (xi , yi , zi , t)]
P
∂aµd
(x, y, z, t) = Fa (x, y, z, t) + µd Ga i=1 N
∂t P
Wa (i, x, y, z, t))
i=1
where a is the quantity being nudged (which can be water vapor mixing ratio, potential temper-
ature, u wind component, or v wind component), µd is the mass of the dry air in the column,
Fa represents the physical tendency terms for a, Ga is the nudging strength for a, N is the
total number of observations, i is the index to the current observation, Wa is the spatiotemporal
weighting function based on the temporal and spatial separation between the observation (lo-
cated at (xi , yi , zi , ti )) and the current model location (x, y, z, t), ao is the observed value of a,
and am (xi , yi , zi , t) is the model value of a interpolated to the observation location. The innova-
tion (ao − am ) for a given observation is calculated once at each time the nudging tendency term
is updated and that innovation is then spread spatially. The nudging strength Ga is the inverse
of the e-folding time of the innovation (assuming that the physical tendency terms [Fa ] were
all zero and thus the observation nudging was the only factor modifying the variable a). Note
that while ARW reads in temperature observations it converts them to potential temperature
for application in observation nudging. Also, surface temperature and surface wind observations
are adjusted to the lowest prognostic level for calculation of the innovation.
103
The user specifies the time period within the model integration over which observation
nudging will be applied. The user also specifies whether the strength at which the nudging is
applied ramps down linearly with time at the end of the observation nudging period. Within the
time period the user specifies for application of observation nudging, individual observations are
applied over a time window of user-specified length centered on the time of the observation. Full
weighting is applied over the half of the window centered at the observation time, and weighting
ramps linearly with time in the quarters of the time window on either end of the time window.
A different time-length can be specified for surface observations than other observations.
Vertically, innovations from surface observations are spread based on user settings and can
be dependent on the depth of the PBL and the current PBL regime. Innovations from multi-
level observations are vertically interpolated in pressure to model levels. For single-level above-
surface observations the innovation is applied from 75 hPa below the observation to 75 hPa
above the observation with weights linearly decreasing with the difference in pressure between
the observation and the level at which the innovation is being applied. Single-level above-surface
observations above (below) the PBL will not be spread below (above) the PBL since errors within
the PBL may not be well-correlated with errors above the PBL. The user can also choose to
prevent any nudging of any observations from being applied within the PBL.
Horizontally, innovations from surface observations are spread along the surface, whereas
innovations from above-surface observations are spread in pressure space (the vertical spreading
described above is applied wherever the innovations are horizontally spread). The user specifies
the horizontal radius of influence for above-surface observations and specifies a factor that this
will be multiplied by to determine the radius of influence for surface observations. The weighting
of the innovation decreases non-linearly with distance from the observation within the radius
of influence. The horizontal radius of influence for above-surface observations increases with
decreasing pressure up until reaching double the specified radius of influence at 500 hPa. An
innovation from a single-level above-surface observation above (below) the PBL will not be
spread horizontally to locations where it would be below (above) the PBL since errors within
the PBL may not be well-correlated with errors above the PBL. Since innovations calculated
in a valley may not be representative of the error on a mountain top, the weighting for surface
observations is decreased based on the difference between the surface pressure at the location of
the observation and at the location the innovation is being applied.
More details regarding observation nudging in ARW can be found in a brief guide on this
topic (http://www2.mmm.ucar.edu/wrf/users/docs/ObsNudgingGuide.pdf).
104
Spectral nudging uses the same input data as for grid-nudging, and nudges model variables
u, v, temperature, and geopotential height. Water mixing ratio nudging was added in Version
4.0. Applications described above for grid-nudging can be considered using spectral nudging
too.
Options to control the grid-nudging are also used for spectral nudging (such as nudging end
time and ramping, nudging strength, controls of nudging in the boundary layer or lower levels
and so on). In addition, there is an option to control the ramping in the vertical between the
layer without nudging to the layer where the nudging is on. The number of waves to nudge in
the x and y directions are user-defined.
105
106
Chapter 11
Data Assimilation
An introduction to the basic ideas of variational data assimilation and the WRF Data Assimila-
tion (WRFDA) system (Barker et al., 2012) is given in this chapter, followed by a brief overview
of recent major improvements to WRFDA. This overview supplements the original description
of the three-dimensional variational (3D-Var) algorithm found in Barker et al. (2004) with sev-
eral important additions and modifications, including a utility gen be for calculating background
error covariances, several new data assimilation algorithms, and new and improved capabilities
to assimilate satellite radiance and radar observations.
11.1 Introduction
The basic goal of any variational data assimilation system is to produce an optimal estimate
of the true atmospheric state at analysis time through iterative solution of a prescribed cost-
function (Ide et al., 1997):
1 1
J(x) = Jb (x) + Jo (x) = (x − xb )T B−1 (x − xb ) + (y − yo )T (E + F)−1 (y − yo ). (11.1)
2 2
The variational problem can be summarized as the iterative minimization of (11.1) to find
the analysis state x that minimizes J(x). This solution represents the a posteriori maximum
likelihood (minimum variance) estimate of the true state of the atmosphere given the two sources
of a priori data: the first guess (or background) xb and observations yo (Lorenc, 1986). The
fit to individual data points is weighted by estimates of their errors: B, E, and F are the back-
ground, observation (instrumental), and representative error covariance matrices, respectively.
The representative error is an estimate of inaccuracies introduced via the observation operator H
used to transform the gridded analysis x to observation space y = H(x) for comparison against
observations. This error will be resolution dependent and may also include a contribution from
approximations (e.g., linearizations) in H.
As described in Barker et al. (2004), the particular variational data assimilation algorithm
adopted in WRFDA is a model-space, incremental formulation of the variational problem. In
this approach, observations, previous forecasts, their errors, and physical laws are combined to
0
produce analysis increments xa , which are added to the first guess xb to provide an updated
analysis.
107
Figure 11.1: Sketch showing the relationship between datasets (circles), and algorithms (rect-
angles) of ARW system.
Figure 11.1 illustrates the relationship between WRFDA, the various datasets, and the other
components of a typical NWP system (here ARW). The WRFDA assimilation proceeds as
described in Barker et al. (2004). A number of recent upgrades to the WRFDA algorithm will
be described in Section 11.2.
The three inputs to WRFDA are:
c) Background error covariances B— used to define the spatial and multivariate response of
108
the analysis to an observation. In variational systems, these covariances are typically calculated
off-line, and significant tuning is required to optimize performance for a particular application
(e.g., Ingleby (2001); Wu et al. (2002)). The amount of work required to do this satisfactorily
is significant, and should not be underestimated. In order to assist the user, the WRFDA
developers supply the following: i) a default set of statistics used for the initial set up of a
domain; ii) a utility gen be (described in Section 11.3) to process ensembles of forecasts into
the appropriate control variable space; and iii) diagnostic routines to assess the accuracy of
observation and background error statistics. These routines include both innovation vector-based
approaches (Hollingsworth and Lonnberg, 1986) and variational tuning approaches (Desroziers
and Ivanov, 2001).
Following assimilation of all data, an analysis xa is produced that must be merged with the
existing lateral boundary conditions xlbc in the UPDATE BC utility (Barker et al. (2003)). At
this stage, the wrfbdy lateral boundary condition files (xlbc ) output of WPS/real is updated to
make the lateral boundaries consistent with the analysis, and surface fields (e.g. SST) are also
updated in the wrfinput analysis file.
109
duced in WRFDA (Zhang et al., 2015). FSO allows the calculation of the forecast impact of
different observations using the adjoint of ARW and WRFDA. Notice that the Lanczos mini-
mization algorithm instead of the Conjugate Gradient Method needs to be used when calculating
FSO with WRFDA.
110
Satellite radiance measurements and RTMs are prone to systematic errors (i.e., biases) that
must be corrected before radiances can be assimilated. Biases typically vary with platform,
instrument, channel, scan angle, and atmospheric conditions. WRFDA adopts the so-called
variational bias correction (VarBC) algorithm for adaptive and online radiance bias correction,
which updates the bias correction coefficients within the linear regression as a part of the vari-
ational minimization (Dee, 2004; Auligne et al., 2007). It is worth mentioning that an offline
VarBC functionality is provided within WRFDA to generate statistics for bias correction co-
efficients without running an actual WRFDA analysis, which is useful when conducting data
assimilation experiments for a short period (Liu et al., 2012).
For radiance assimilation under clear-sky condition, two advanced cloud detection methods
for hyperspectral infrared sensors such as AIRS and IASI were implemented in WRFDA (Xu et
al., 2013, 2014, 2015; Auligne, 2014a,b), which allows the assimilation of the channels peaking
above the cloud. In version 3.9, all-sky radiance assimilation capability was introduced for
AMSR2 radiance data (Yang et al., 2016) with the development of the so-called “symmetric
error model” for all-sky data (Geer and Bauer, 2011).
x0 = Uv = Up Uv Uh v. (11.2)
The expansion U = Up Uv Uh represents the various stages of covariance modeling: horizontal
correlations Uh through recursive filter, vertical covariances Uv through EOF decomposation,
and multivariate covariances Up through statistical regression.
The components of v are chosen so that their error cross-correlations are negligible, thus
permitting the matrix B to be block-diagonalized. There are three choices of control variables
in WRFDA through a namelist parameter “cv options”:
111
• cv options=5: Streamfunction ψ 0 , unbalanced velocity potential χ0u , unbalanced tempera-
ture Tu0 , pseudo relative humidity rh0 , unbalanced surface pressure ps0u . This option allows
the statistical cross-correlation between mass and wind field and univariate rh analysis
(i.e., no cross correlation between rh and other variables).
• cv options=6 (Chen et al., 2013): Streamfunction ψ 0 , unbalanced velocity potential χ0u ,
unbalanced temperature Tu0 , unbalanced pseudo relative humidity rh0u , unbalanced surface
pressure ps0u . In addition to mass-wind balance, this option allows the cross-correlation
between moisture and other variables, which may be more important in tropical regions.
• cv options=7 (Sun et al., 2015): zonal wind u0 , meridional wind v 0 , temperature T 0 , pseudo
relative humidity rh0 , surface pressure ps0 . There is no cross-correlation between those
variables (so completely univariate analysis).
a) Introduced in Version 3.8, the “weak penalty constraint” (WPEC) option (Li et al., 2015)
aims to enforce quasi-gradient balance on the analysis. It was designed with the specific purpose
of improving the assimilation of radar data within tropical cyclones, but may be useful for other
weather phenomena of similar scales. It can be used with 3D-Var or hybrid-3DEnVar (4D-Var
is not compatible with this capability).
b) In Version 4.0, a divergence constraint (DIVC) term was added to model the correlation
between u and v. The DIVC was implemented by adding a term in the cost function that
constrains the horizontal divergence (Tong et al., 2016).
c) Introduced in Version 4.0, the large scale analysis constraint (LSAC) option (Vendrasco et
al., 2016) ensures that the convective-scale analysis does not distort the underlining large-scale
balance and eliminates possible large-scale bias in the ARW background. The global analysis
or forecast, such as that from GFS or FNL, is treated as bogus observations and assimilated via
112
WRFDA. The input data for LSAC is prepared using WPS/REAL, the same as preparing for
ARW input data.
d) A new option to assimilate GPS Radio Occultation (GPSRO) refractivity data using the
GPS Excess Phase (GPSEPH) nonlocal operator (Chen et al., 2009) is added in Version 4.0
with a parallelization strategy described in (Zhang et al., 2014b).
e) Since Version 3.5, WRFDA can assimilate wind observations in terms of wind speed and
direction (Huang et al., 2013; Gao et al., 2015), in addition to originally assimilating u and v
components.
113
error is not known in reality, but is assumed to be statistically well-represented by a model state
perturbation x0 . In the standard NMC-method (Parrish and Derber, 1992), the perturbation
x0 is given by the difference between two forecasts (e.g., 24 hour minus 12 hour) verifying
at the same time. Climatological estimates of background error may then be obtained by
averaging such forecast differences over a period of time (e.g., one month). An alternative
strategy proposed by (Fisher, 2003) makes use of ensemble forecast output, defining the x0
vectors as ensemble perturbations (ensemble minus ensemble mean). In either approach, the
end result is an ensemble of model perturbation vectors from which estimates of background error
may be derived. The gen be utility has been designed to work with either forecast difference,
or ensemble-based, perturbations. Using the NMC-method, x0 = xT2 − xT1 where T 2 and T 1
are the forecast difference times (e.g., 48h minus 24h for global, 24h minus 12h for regional).
Alternatively, for an ensemble-based approach, xk 0 = xk − x̄, where the overbar is an average
over ensemble members k = 1, ne . The total number of binary files produced by this stage is
nf × ne where nf is the number of forecast times used (e.g., for 30 days with forecast every
12 hours, nf = 60). Using the NMC-method, ne = 1 (1 forecast difference per time). For
ensemble-based statistics, ne is the number of ensemble members.
As described above, the WRFDA background error covariances are specified not in model
space x0 , but in a control variable space v, which is related to the model variables (e.g., wind
components, temperature, humidity, and surface pressure) via the control variable transform
defined in (11.2). Both (11.2) and its adjoint are required in WRFDA. To enable this, the
(offline) background error utility is used to compute components of the forecast error covariance
matrix modeled within the U transform. This process is described in the following subsections.
The background error covariance generation code gen be is designed to process data from a
variety of regional/global models (e.g., ARW, MM5, KMA global model, NFS, etc.), and process
it in order to provide error covariance statistics for use in variational data assimilation systems.
The initial, model-dependent “stage 0” is illustrated in Fig. 11.2.
Alternative models use different grids, variables, data formats, etc., and so initial converters
are required to transform model output into a set of standard perturbation fields (and metadata),
and to output them in a standard binary format for further, model-independent processing. The
standard grid-point fields are as follows.
• Perturbations: Streamfunction ψ 0 (i, j, k), velocity potential χ0 (i, j, k), temperature T 0 (i, j, k),
relative humidity r0 (i, j, k), surface pressure p0s (i, j).
• Full-fields: height z(i, j, k), latitude φ(i, j). (These are required for the production of
background error statistics stored in terms of physics variables, rather than tied to a
specific grid. This flexibility is included in gen be through a namelist option to define the
bins over which data is averaged in a variety of ways (e.g., latitude height, grid points).
Land-sea and orographic effects may also be represented in this way.
In general, the stage 0 converters are developed and maintained by those supporting indi-
vidual models. Only the WRF-to-standard-fields converter gen be stage0 wrf is maintained and
supported by ARW effort.
114
Figure 11.2: Sketch of the role of Stage 0 converters in transforming model-specific data (e.g.,
ARW, KMA global model, etc.) to standard perturbation fields and relevant metadata (e.g.,
latitude, height, land/sea, etc.).
The summation over the vertical index k1 relates to the integral (hydrostatic) relationship
between mass fields and the wind field. By default, the regression coefficients c, G, and W do
not vary horizontally, however options exists to relax this assumption via the bin type namelist
variable in order to allow representation of differences between, for example, polar, mid-latitude,
and tropical dynamical and physical processes. The scalar coefficient c used to estimate velocity
potential errors from those of streamfunction is permitted to vary with model level in order to
represent, for example, the impact of boundary-layer physics. Latitudinal/height smoothing of
115
the resulting coefficients may be optionally performed to avoid artificial discontinuities at the
edges of latitude/height boxes.
Having computed regression coefficients, the unbalancedP components of the fields are cal-
culated
P as χ u (k) = χ(k) − c(k)ψ(k), Tu (k) = T (k) − k1 G(k1, k)ψ(k1), and psu = ps −
k1 W (k1)ψ(k1). These fields are output for the subsequent calculation of the spatial covari-
ances as described below.
116
Appendix A
Physical Constants
π = 3.1415926 Pi
k = 0.4 Von Karman constant
re = 6.370 × 106 m Radius of earth
g = 9.81 m s−2 Acceleration due to gravity
Ωe = 7.2921 × 10−5 s−1 Angular rotation rate of the earth
σB = 5.67051 × 10−8 W m−2 K−4 Stefan − Boltzmann constant
Rd = 287 J kg−1 K−1 Gas constant for dry air
Rv = 461.6 J kg−1 K−1 Gas constant for water vapor
cp = 7 × Rd /2 J kg−1 K−1 Specific heat of dry air at constant pressure
cv = cp − Rd J kg−1 K−1 Specific heat of dry air at constant volume
cpv = 4 × Rv J kg−1 K−1 Specific heat of water vapor at constant pressure
cvv = cpv − Rv J kg−1 K−1 Specific heat of water vapor at constant volume
cliq = 4190 J kg−1 K−1 Specific heat capacity of water
cice = 2106 J kg−1 K−1 Specific heat capacity of ice
Lv = 2.5 × 106 J kg−1 Latent heat of vaporization
Ls = 2.85 × 106 J kg−1 Latent heat of sublimation
Lf = 3.50 × 105 J kg−1 Latent heat of fusion
ρw = 1.0 × 103 kg m−3 Density of liquid water
117
118
Appendix B
List of Symbols
Symbols used in this document are listed in alphabatical order in this appendix.
Symbols Definition
a generic variable
A coefficient (Chapter 4), base-state lapse rate constant (Chapter 5)
B Vertical profile in hydrid coordinate definition (Chapter 2)
B background error covariance matrix
c scalar coefficient
cs speed of sound
Ck a constant used in TKE closure
Cr Courant number
Crmax maximum Courant number
Crtheory Courant number from Table 3.1
Crβ activation Courant number in vertical velocity damping
Cs a constant used in eddy viscosity calculation
D deformation
Dnm deformation tensor, where n, m = 1, 2 and 3
e cosine component of the Coriolis term (Chapters 2, 3); turbulent kinetic energy
(Chapter 4)
E observation error covariance matrix
f sine component of the Coriolis term
F forcing terms for U , V , W , Θ and Qm
F representivity error covariance matrix
FXcor Coriolis forcing terms for X = U , V , and W
F1,2 coefficients for weighting functions in specified boundary condition
g acceleration due to gravity
Gk regression coefficient
H observation operator
J cost function
Kdh,dv horizontal and vertical eddy viscosity for gravity wave absorbing layer
Kh,v horizontal and vertical eddy viscosities
119
Symbols Definition
120
Symbols Definition
Wk regression coefficient
z height
zd depth of damping layer
ztop height of model top
α inverse density of air
α0 perturbation inverse density of air
ᾱ inverse density of air for the reference state
αd inverse density of dry air
αr local rotation angle between y-axis and the meridian
β off-centering coefficient for semi-implicit acoustic step
γ ratio of heat capacities for dry air at constant pressure and volume
γd divergence damping coefficient
γe external mode damping coefficient
γg damping coefficient for upper boundary gravity wave absorbing layer
γr Rayleigh damping coefficient
molecular weight of water over the molecular weight of dry air (Chapter 4);
true background error (Chapter 11)
η terrain-following hybrid sigma-pressure vertical coordinate
η̇ contravariant ‘vertical’ velocity or coordinate velocity
θ potential temperature
θe equivalent potential temperature
θm moist potential temperature
Θm coupled moist potential temperature
µ̄d reference state vertical coordinate metric
µd vertical coordinate metric
τ acoustic time (Chapter 3), vertical structure function for Rayleigh damping
(Chapter 4)
τnm stress tensor (Chapter 4) where n.m = 1, 2 and 3
∆τ acoustic time step
φ geopotential (Chapters 2, 3, 5); latitude (Chapter 11)
φ̄ geopotential for reference state
φ0 perturbation geopotential
Φ generic prognostic variable (coupled)
ψ generic variable (Chapter 6)
ψ0 streamfunction increment
χ0 velocity potential increment
ω same as η̇
Ω coupled coordinate velocity
Ωe angular rotation rate of the earth
121
Subscripts/Superscripts Definition
()d dry
()h hydrostatic
()0 base state sea-level constant
() reference state
()0 perturbation from reference state
∗
()t full value at a Runge-Kutta step
()00 perturbation from Runge-Kutta step value in acoustic steps
122
123
Appendix C
Acronyms
125
126
Bibliography
Auligne, T. D., 2014: Multivariate minimum residual method for cloud retrieval. Part I: Theo-
retical aspects and simulated observation experiments. Mon. Wea. Rev., 142, 4383–4398.
Auligne, T. D., 2014: Multivariate minimum residual method for cloud retrieval. Part II: Real
observations experiments. Mon. Wea. Rev., 142, 4399–4415.
Auligne, T., A. P. McNally, and D. P. Dee, 2007: Adaptive bias correction for satellite data in
a numerical weather prediction system. Q. J. R. Meteorol. Soc., 133, 631–642.
Bae, S. Y., S.-Y. Hong, and W.-K. Tao, 2018: Development of a single-moment cloud micro-
physics scheme with prognostic hail for the Weather Research and Forecasting (WRF) model.
Asia-Pac. J. Atmos. Sci., https://doi.org/10.1007/s13143-018-0066-3.
Baek, S., 2017: A revised radiation package of G-packed McICA and two-stream approximation:
Performance evaluation in a global weather forecasting model. J. Adv. Model. Earth Syst., 9,
1–13, doi:10.1002/2017MS000994.
Ban, J., Z. Liu, X. Zhang, X.-Y. Huang, and H. Wang, 2017: Precipitation data assimila-
tion in WRFDA 4D-Var: implementation and application to convection-permitting forecasts
over United States. Tellus A: Dynamic Meteorology and Oceanography, 69:1, 1368310, DOI:
10.1080/16000870.2017.1368310.
Barker, D. M., W. Huang, Y.-R. Guo, and A. Bourgeois, 2003: A Description of the Ad-
vanced Research WRF Version 3. NCAR Tech Note, NCAR/TN-475+STR, pp. [Available
from UCAR Communications, P.O. Box 3000, Boulder, CO, 80307.]
127
Barker, D. M., W. Huang, Y.-R. Guo, A. Bourgeois, and X. N. Xiao, 2004: A Three-Dimensional
Variational Data Assimilation System for MM5: Implementation and Initial Results Mon.
Wea. Rev., 132, 897–914.
Beljaars, A.C.M., 1994: The parameterization of surface fluxes in large-scale models under free
convection, Quart. J. Roy. Meteor. Soc., 121, 255–270.
Berg, L.K., W.I. Gustafson, E.I. Kassianov, E. I., and L. Deng (2013), Evaluation of a modified
scheme for shallow convection: Implementation of CuP and case studies, Mon. Wea. Rev.,
141, 134-147, 10.1175/mwr-d-12-00136.1.
Berg, L.K, M. Shrivastava, R.C. Easter, J.D. Fast, E.G. Chapman, and Y. Liu, 2015: A new
WRF-Chem treatment for studying regional scale impacts of cloud-aerosol interactions in
parameterized cumuli. Geophys. Model Devel., 8, 409-429. doi:10.5194/gmd-8-409-2015.
Berner, J., K. R. Fossell, S.-Y. Ha, J. P. Hacker, and C. Snyder, 2015: Increasing the skill of
probabilistic forecasts: Understanding performance improvements from model-error represen-
tations. Mon. Wea. Rev., 143, 1295–1320.
Berner, J., S.-Y. Ha, J. P. Hacker, A. Fournier, and C. Snyder, 2011: Model uncertainty in a
mesoscale ensemble prediction system: Stochastic versus multi-physics representations. Mon.
Wea. Rev., 139, 1972–1995.
Berner, J., G. Shutts, M. Leutbecher, and T. Palmer, 2009: A spectral stochastic kinetic energy
backscatter scheme and its impact on flow-dependent predictability in the ECMWF ensemble
prediction system. J. Atmos. Sci., 66, 603–626.
Berner, J., and Coauthors, 2017: Stochastic parameterization: Toward a new view of weather
and climate models. Bulletin of the American Meteorological Society, 98 (3), 565–588.
Betts, A. K., 1986: A new convective adjustment scheme. Part I: Observational and theoretical
basis. Quart. J. Roy. Meteor. Soc., 112, 677–691.
Betts, A. K., and M. J. Miller, 1986: A new convective adjustment scheme. Part II: Single column
tests using GATE wave, BOMEX, and arctic air-mass data sets. Quart. J. Roy. Meteor. Soc.,
112, 693–709.
Bowler, N. E., A. Arribas, K. R. Mylne, K. B. Robertson, and S. E. Beare, 2008: The MOGREPS
short-range ensemble prediction system. Quart. J. Roy. Meteor. Soc., 134, 703–722.
128
Bretherton, C. S., and S. Park, 2009: A new moist turbulence parameterization in the Commu-
nity Atmosphere Model. J. Climate, 22, 3422–3448.
Buizza, R., M. Miller, and T. N. Palmer, 1999: Stochastic representation of model uncertainties
in the ECMWF Ensemble Prediction System.
Cavallo, S. M., J. Dudhia and C. Snyder, 2011: A multi-layer upper boundary condition for
longwave radiative flux to correct temperature biases in a mesoscale model. Mon. Wea. Rev.,
139, 1952–1959.
Chen, F., and J. Dudhia, 2001: Coupling an advanced land-surface/ hydrology model with the
Penn State/ NCAR MM5 modeling system. Part I: Model description and implementation.
Mon. Wea. Rev., 129, 569–585.
Chen, S. Y., C. Y. Huang, Y. H. Kuo, Y. R. Guo, and S. Sokolovskiy, 2009: Assimilation of GPS
refractivity from FORMOSAT-3/COSMIC using a nonlocal operator with WRF 3DVAR and
its impact on the prediction of a typhoon event. Terr. Atmos. Ocean. Sci., 20, 133–154.
Chen, F., M. Tewari, H. Kusaka, and T. T. Warner, 2006: Current status of urban modeling in
the community Weather Research and Forecast (WRF) model. Joint with Sixth Symposium
on the Urban Environment and AMS Forum: Managing our Physical and Natural Resources:
Successes and Challenges, Atlanta, GA, USA, Amer. Meteor. Soc., CD-ROM. J1.4.
Chen, S.-H., and W.-Y. Sun, 2002: A one-dimensional time dependent cloud model. J. Meteor.
Soc. Japan, 80, 99–118.
Chen, Y., S. R. H. Rizvi, X.-Y. Huang, J. Min, and X. Zhang, 2013: Balance characteristics
of multivariate background error covariances and their impact on analyses and forecasts in
tropical and Arctic regions. Meteor. Atmos. Phys., 121, 79–98.
Choi H., and S. Hong, 2015: An updated subgrid orographic parameterization for global atmo-
spheric forecast models. J. Geophys. Res., 120, 12445–12457.
Chou M.-D., and M. J. Suarez, 1994: An efficient thermal infrared radiation parameterization
for use in general circulation models. NASA Tech. Memo. 104606, 3, 85pp.
Chou M.-D., and M. J. Suarez, 1999: A solar radiation parameterization for atmospheric studies.
NASA Tech. Memo. 104606, 15, 40pp.
Chou, M. D., M. J. Suarez, X. Z. Liang, and M. M. H. Yan, 2001: A thermal infrared radiation
parameterization for atmospheric studies. NASA Tech. Memo. 104606, 19, 68pp.
Collins, W.D. et al., 2004: Description of the NCAR Community Atmosphere Model (CAM
3.0), NCAR Technical Note, NCAR/TN-464+STR, 226pp.
129
Cooper, W. A., 1986: Ice initiation in natural clouds. Precipitation Enhancement — A Scientific
Challenge, Meteor. Monogr., No. 43, Amer. Met. Soc., 29–32.
Daley, R., and E. Barker, 2001: NAVDAS: Formulation and Diagnostics. Mon. Wea. Rev., 129,
869–883.
Daniels, M., K. Lundquist, J. Mirocha, D. Wiersema, and F. Chow, 2016. A new vertical grid
nesting capability in the Weather Research and Forecasting (WRF) model. Mon. Wea. Rev.,
144, 3725—3747. DOI:10.1175/MWR-D-16-0049.1
Davies, H. C., and R. E. Turner, 1977: Updating prediction models by dynamical relaxation:
An examination of the technique. Quart. J. Roy. Meteor. Soc., 103, 225–245.
Deardorff, J. W., 1972 Numerical investigation of neutral and unstable planetary boundary
layers, J. Atmos. Sci., 29, 91–115.
Dee, D. P.,2004: Variational bias correction of radiance data in the ECMWF system. em In:
Proc. of the ECMWF Workshop on Assimilation of High Spectral Resolution Sounders in
NWP, Reading, U. K., 28 June–1 July.
Desroziers, G., 1997: A coordinate change for data assimilation in spherical geometry of frontal
structures. Mon. Wea. Rev., 125, 3030–3039.
Desroziers, G., and S. Ivanov, 2001: Diagnosis and adaptive tuning of observation-error param-
eters in a variational assimilation. Quart. J. Roy. Meteor. Soc., 127, 1433–1452.
Dudhia, J., 1989: Numerical study of convection observed during the winter monsoon experiment
using a mesoscale two-dimensional model, J. Atmos. Sci., 46, 3077–3107.
Dudhia, J., S.-Y. Hong, and K.-S. Lim, 2008: A new method for representing mixed-phase
particle fall speeds in bulk microphysics parameterizations. J. Met. Soc. Japan, 86, 33–44.
Durran, R. D., and J. B. Klemp, 1983: A compressible model for the simulation of moist
mountain waves, Mon. Wea. Rev., 111, 2341–2361.
Dyer, A. J., and B. B. Hicks, 1970: Flux-gradient relationships in the constant flux layer, Quart.
J. Roy. Meteor. Soc., 96, 715–721.
130
Easter, R. C., 1993: Two modified versions of Bott’s positive definite numerical advection
scheme. Mon. Wea. Rev., 121, 297–304.
Fisher, M., 2003: Background error covariance modeling. Seminar on Recent Development in
Data Assimilation for Atmosphere and Ocean, 45–63, ECMWF.
Fletcher, N. H., 1962: The Physics Of Rain Clouds. Cambridge University Press, 386 pp.
Fu, Q., and K. N. Liou, 1992: On the correlated k-distribution method for radiative transfer in
nonhomogeneous atmospheres. J. Atmos. Sci., 49, 2139–2156.
Gao, F., X.-Y. Huang, N. A. Jacobs, and H. Wang, 2015: Assimilation of wind speed and
direction observations: results from real observation experiments. Tellus A, 67, 27132,
doi:10.3402/tellusa.v67.27132.
Geer, A. J. and Bauer, P., 2011: Observation errors in all-sky data assimilation. Q. J. Roy.
Meteorol. Soc., 137, 2024–2037.
Glotfelty, T., K. Alapaty, J. He, P. Hawbecker, X. Song, and G. Zhang, 2019: The Weather
Research and Forecasting Model with Aerosol Cloud Interactions (WRF-ACI): Development,
Evaluation, and Initial Application. Mon. Wea. Rev., in press.
Grell, G. A., and D. Devenyi, 2002: A generalized approach to parameterizing convection com-
bining ensemble and data assimilation techniques. Geophys. Res. Lett., 29(14), Article 1693.
Grell, G.A., S.E. Peckham, R. Schmitz, S.A. McKeen, G. Frost, W.C. Skamarock and B. Eder,
2005: Fully coupled online chemistry within the WRF model. Atmos. Environ., 39, 6957-6975.
Grell, G. A. and Freitas, S. R., 2014: A scale and aerosol aware stochastic convective pa-
rameterization for weather and air quality modeling. Atmos. Chem. Phys., 14, 5233-5250,
doi:10.5194/acp-14-5233-2014.
Grenier, H., and C. S. Bretherton, 2001: A moist PBL parameterization for large-scale models
and its application to subtropical cloud-topped marine boundary layers. Mon. Wea. Rev., 129,
357–377.
131
Gu, Y., K. N. Liou, S. C. Ou, and R. Fovell, 2011: Cirrus cloud simulations using WRF with
improved radiation parameterization and increased vertical resolution. J. Geophys. Res., 116,
D06119.
Hacker, J. P., C. Snyder, S.-Y. Ha, and M. Pocernich, 2011: Linear and nonlinear response to
parameter variations in a mesoscale model. Tellus A, 63, 429–444.
Haltiner, G. J., and R. T. Williams, 1980: Numerical prediction and dynamic meteorology. John
Wiley & Sons, Inc., 477pp.
Han, Jongil and HuaLu Pan, 2011: Revision of convection and vertical diffusion schemes in
the NCEP Global Forecast System. Wea. Forecasting, 26, 520533. doi:10.1175/WAF-D-10-
05038.1.
Hollingsworth, A., and P. Lonnberg, 1986: The Statistical Structure Of Short-Range Forecast
Errors As Determined From Radiosonde Data. Part I: The Wind Field. Tellus, 38A, 111–136.
Holtslag, A. A. M. and B. A. Boville, 1993: Local versus non-local boundary layer diffusion in
a global climate model, J. Climate, 6, 1825–1842.
Hong, S.-Y., 2007: Stable Boundary Layer Mixing in a Vertical Diffusion Scheme. The Korea
Meteor. Soc., Fall conference, Seoul, Korea, Oct. 25-26.
Hong, S.-Y. and J.-O. J. Lim, 2006: The WRF Single-Moment 6-Class Microphysics Scheme
(WSM6), J. Korean Meteor. Soc., 42, 129–151.
Hong, S.-Y. and H.-L. Pan, 1996: Nonlocal boundary layer vertical diffusion in a medium-range
forecast model, Mon. Wea. Rev., 124, 2322–2339.
Hong, S.-Y., H.-M. H. Juang, and Q. Zhao, 1998: Implementation of prognostic cloud scheme
for a regional spectral model, Mon. Wea. Rev., 126, 2621–2639.
Hong, S.-Y., J. Dudhia, and S.-H. Chen, 2004: A Revised Approach to Ice Microphysical Pro-
cesses for the Bulk Parameterization of Clouds and Precipitation, Mon. Wea. Rev., 132,
103–120.
Hong, S.-Y., Y. Noh, and J. Dudhia, 2006: A new vertical diffusion package with an explicit
treatment of entrainment processes. Mon. Wea. Rev., 134, 2318–2341.
Hong, S.-Y. and J.-H. Jang, 2018: Impacts of shallow convection processes on a simulated
Boreal summer climatology in a global atmospheric model. Asia-Pacific J. Atmos. Sci. (2018)
54(Suppl 1): 361. doi:10.1007/s13143-018-0013-3.
Huang, X.-Y. and P. Lynch, 1993: Diabatic Digital-Filtering Initialization: Application to the
HIRLAM Model., Mon. Wea. Rev., 121, 589–603.
132
Huang, X.-Y., F. Gao, N. A. Jacobs, and H. Wang, 2013: Assimilation of wind speed and
direction observations: a new formulation and results from idealised experiments. Tellus A,
65, 19936, doi:10.3402/tellusa.v65i0.19936.
Ide, K., P. Courtier, M. Ghil, and A. C. Lorenc, 1997: Unified notation for data assimilation:
Operational, sequential and variational. J. Met. Soc. Japan, 75, 181–189.
Ingleby, N. B., 2001: The statistical structure of forecast errors and its representation in the
Met. Office global 3-D variational data assimilation scheme. Quart. J. Roy. Meteor. Soc., 127,
209–232.
Jankov, I., and Coauthors, 2017: A performance comparison between multiphysics and stochas-
tic approaches within a north american rap ensemble. Monthly Weather Review, 145, 1161–
1179.
Janjic, Z. I., 1990: The step-mountain coordinate: physical package, Mon. Wea. Rev., 118,
1429–1443.
Janjic, Z. I., 1994: The step-mountain eta coordinate model: further developments of the con-
vection, viscous sublayer and turbulence closure schemes, Mon. Wea. Rev., 122, 927–945.
Janjic, Z. I., 1996: The surface layer in the NCEP Eta Model, Eleventh Conference on Numerical
Weather Prediction, Norfolk, VA, 19–23 August; Amer. Meteor. Soc., Boston, MA, 354–355.
Janjic, Z. I., 2000: Comments on ”Development and Evaluation of a Convection Scheme for Use
in Climate Models”, J. Atmos. Sci., 57, p. 3686.
Janjic, Z. I., 2002: Nonsingular Implementation of the Mellor–Yamada Level 2.5 Scheme in the
NCEP Meso model, NCEP Office Note, No. 437, 61 pp.
Jimenez, P. A., and J. Dudhia, 2012: Improving the representation of resolved and unresolved
topographic effects on surface wind in the WRF model. J. Appl. Meteor. Climatol., 51, 300–
316.
Kain, J. S., and J. M. Fritsch, 1990: A one-dimensional entraining/ detraining plume model
and its application in convective parameterization, J. Atmos. Sci., 47, 2784–2802.
133
Kain, J. S., and J. M. Fritsch, 1993: Convective parameterization for mesoscale models: The
Kain-Fritcsh scheme, The representation of cumulus convection in numerical models, K. A.
Emanuel and D.J. Raymond, Eds., Amer. Meteor. Soc., 246 pp.
Kain, J. S., 2004: The Kain-Fritsch convective parameterization: An update. J. Appl. Meteor.,
43, 170–181.
Kessler, E., 1969: On the distribution and continuity of water substance in atmospheric circu-
lation, Meteor. Monogr., 32, Amer. Meteor. Soc., 84 pp.
Khain, A., B. Lynn, and J. Dudhia, 2010: Aerosol effects on intensity of landfalling hurricanes
as seen from simulations with the WRF model with spectral bin microphysics. J. Atmos. Sci.,
67, 365–384.
Khain, A., A. Pokrovsky, M. Pinsky, A. Seifert, and V. Phillips, 2004: Simulation of effects
of atmospheric aerosols on deep turbulent convective clouds using a spectral microphysics
mixed-phase cumulus cloud model. Part I: model description and possible applications. J.
Atmos. Sci., 61, 2963–2982.
Klemp, J. B., and D. K. Lilly, 1978: Numerical simulation of hydrostatic mountain waves, J.
Atmos. Sci., 35, 78–107.
Klemp, J. B., and R. Wilhelmson, 1978: The simulation of three-dimensional convective storm
dynamics, J. Atmos. Sci., 35, 1070–1096.
Klemp, J. B., W. C. Skamarock, and J. Dudhia, 2007: Conservative split-explicit time in-
tegration methods for the compressible nonhydrostatic equations, Mon. Wea. Rev.., 135,
2897–2913.
Klemp, J. B., J. Dudhia, and A. Hassiotis, 2008: An Upper Gravity Wave Absorbing Layer for
NWP Applications. Mon. Wea. Rev.., 136, 3987–4004.
Knievel, J. C., G. H. Bryan, and J. P. Hacker, 2007: Explicit numerical diffusion in the WRF
Model. Mon. Wea. Rev., 135, 3808–3824.
Kusaka, H. and F. Kimura, 2004: Coupling a single-layer urban canopy model with a simple
atmospheric model: Impact on urban heat island simulation for an idealized case. J. Meteor.
Soc. Japan, 82, 67–80.
Kusaka, H., H. Kondo, Y. Kikegawa, and F. Kimura, 2001: A simple single-layer urban canopy
model for atmospheric models: Comparison with multi-layer and slab models. Bound.-Layer
Meteor., 101, 329–358.
Kwon, Y.-C. and S.-Y. Hong, 2017: A mass-flux cumulus parameterization scheme across gray-
zone resolutoins. Mon. Wea. Rev. 145, 583-598, doi:0.1175/MWR-D-16-0034.1
Lacis, A. A., and J. E. Hansen, 1974: A parameterization for the absorption of solar radiation
in the earth’s atmosphere. J. Atmos. Sci., 31, 118–133.
134
Lang, S., W.-K. Tao, J.-D. Chern, D. Wu, and X. Li, 2014: Benefits of a 4th ice class in
the simulated radar reflectivities of convective systems using a bulk microphysics scheme. J.
Atmos. Sci., 71, 3583–3612, https://doi.org/10.1175/JAS-D-13-0330.1
Laprise R., 1992: The Euler Equations of motion with hydrostatic pressure as as independent
variable, Mon. Wea. Rev., 120, 197–207.
Lawrence, D. M., et al., 2011: Parameterization improvements and functional and structural
advances in Version 4 of the Community Land Model. J. Adv. Model. Earth Syst., 3, M03001.
Lee, C.-Y. and S. S. Chen, 2012: Symmetric and asymmetric structures of hurricane boundary
layer in coupled atmosphere-wave-ocean models and observations. J. Atmos. Sci., 69, 3576–
3594.
Lee, M.-S., D. Barker, W. Huang and Y.-H. Kuo, 2004: First Guess at Appropriate Time
(FGAT) with WRF 3DVAR. WRF/MM5 Users Workshop, 22–25 June 2004, Boulder, Col-
orado.
Li, D., E. Bou-Zeid, M. Barlage, F. Chen, and J. A. Smith, 2013: Development and Evaluation
of a Mosaic Approach in the WRF-Noah Framework. J. Geophys. Res., 118, 11918–11935.
Li, X., J. Ming, M. Xue, Y. Wang, and K. Zhao, 2015: Implementation of a dynamic equation
constraint based on the steady state momentum equations within the WRF hybrid ensemble-
3DVar data assimilation system and test with radar T-TREC wind assimilation for tropical
Cyclone Chanthu (2010). J. Geophys. Res. Atmos., 120, 4017–4039.
Lim, K.-S. S., and S.-Y. Hong, 2010: Development of an effective double-moment cloud mi-
crophysics scheme with prognostic cloud condensation nuclei (CCN) for weather and climate
models. Mon. Wea. Rev., 138, 1587–1612.
Lin, Y., and B. A. Colle, 2011: A new bulk microphysical scheme that includes riming intensity
and temperature-dependent ice characteristics. Mon. Wea. Rev. , 139, 1013–1035.
Lin, Y.-L., R. D. Farley, and H. D. Orville, 1983: Bulk parameterization of the snow field in a
cloud model. J. Climate Appl. Meteor., 22, 1065–1092.
135
Liu Z. and Barker D. M., 2006: Radiance Assimilation in WRF-Var: Implementation and Initial
Results, Preprints, 7th WRF Users’ Workshop, Boulder, CO, National Center for Atmospheric
Research.
Liu, Z., C. S. Schwartz, C. Snyder, and S.-Y. Ha, 2012: Impact of assimilating AMSU-A radi-
ances on forecasts of 2008 Atlantic tropical cyclones initialized with a limited-area ensemble
Kalman filter. Mon. Wea. Rev., 140, 4017–4034.
Lorenc, A. C., 1986: Analysis methods for numerical weather prediction. Quart. J. Roy. Meteor.
Soc., 112, 1177–1194.
Lorente-Plazas, R., P. A. Jimenez, J. Dudhia, and J. P. Montavez, 2016: Evaluating and im-
proving the impact of the atmospheric stability and orography on surface winds in the WRF
model. Mon. Wea. Rev., 144, 2685–2693.
Lynch, P., and X.-Y. Huang, 1992: Initialization of the HIRLAM Model Using a Digital Filter.
Mon. Wea. Rev., 120, 1019–1034.
Lynch, P., and X.-Y. Huang, 1994: Diabatic initialization using recursive filters. Tellus, 46A,
583–597.
Lynch, P., 1997: The Dolph-Chebyshev Window: A Simple Optimal Filter. Mon. Wea. Rev.,
125, 655–660.
Ma, L.-M. and Z.-M. Tan, 2009: Improving the behavior of the cumulus parameterization
for tropical cyclone prediction: Convection trigger. Atmospheric Research, 92, 190–211,
doi:10.1016/j.atmosres.2008.09.022.
Mahalov, A., and M. Moustaoui, 2009: Vertically nested nonhydrostatic model for multiscale
resolution of flows in the upper troposphere and lower stratosphere. J. Comput. Phys., 228,
1294–1311, doi:10.1016/j.jcp.2008.10.030.
Martilli A, Clappier A, and Rotach M.W., 2002: An urban surface exchange parameterization
for mesoscale models. Bound.-Layer Meteorol., 104, 261–304.
Matsui, T., S. Q. Zhang, W.-K. Tao, S. Lang, C. Ichoku, and C. Peters-Lidard, 2018: Im-
pact of Radiation Frequency, Precipitation Radiative Forcing, and Radiation Column Ag-
gregation on Convection-Permitting West African Monsoon Simulations. Clim. Dyn., 1-21,
https://doi.org/10.1007/s00382-018-4187-2
McCumber, M., W.-K. Tao, J. Simpson, R. Penc, and S.-T. Soong, 1991: Comparison of ice-
phase microphysical parameterization schemes using numerical simulations of tropical con-
vection. J. Appl. Meteor., 30, 985–1004.
Mellor, G. L., and T. Yamada, 1982: Development of a turbulence closure model for geophysical
fluid problems. Rev. Geophys. Space Phys., 20, 851–875.
136
Michalakes, J., J. Dudhia, D. Gill, J. Klemp, and W. Skamarock, 1999: Design of a next-
generation regional weather research and forecast model, Towards Teracomputing, World Sci-
entific, River Edge, New Jersey, 117–124.
Miguez-Macho, G., G. L. Stenchikov, and A. Robock, 2004: Spectral nudging to eliminate the
effects of domain position and geometry in regional climate model simulations. J. Geophys.
Res., 109, D13104. doi:10.1029/2003JD004495.
Monin, A.S. and A.M. Obukhov, 1954: Basic laws of turbulent mixing in the surface layer of
the atmosphere. Contrib. Geophys. Inst. Acad. Sci., USSR, (151), 163–187 (in Russian).
Morrison, H., and J. O. Pinto, 2006: Intercomparison of bulk microphysics schemes in mesoscale
simulations of springtime Arctic mixed-phase stratiform clouds. Mon. Wea. Rev., 134, 1880–
1900.
Morrison, H., and A. Gettelman, 2008: A New Two-Moment Bulk Stratiform Cloud Micro-
physics Scheme in the Community Atmosphere Model, Version 3 (CAM3). Part I: Description
and Numerical Tests. J. Clim.. 21, 3642–3659.
Morrison, H. and J.A. Milbrandt, 2015: Parameterization of ice microphysics based on the
prediction of bulk particle properties. Part 1: Scheme description and idealized tests. J.
Atmos. Sci., 72, 287–311.
Morrison, H., G. Thompson, and V. Tatarskii, 2008: Impact of cloud microphysics on the
development of trailing stratiform precipitation in a simulated squall line: Comparison of
one- and two-moment schemes. Mon. Wea. Rev. , 137, 991–1007.
Nakanishi, M., and H. Niino, 2006: An improved Mellor-Yamada level 3 model: its numerical
stability and application to a regional prediction of advecting fog. Bound. Layer Meteor., 119,
397–407.
Nakanishi, M., and H. Niino, 2009: Development of an improved turbulence closure model for
the atmospheric boundary layer. J. Meteor. Soc. Japan, 87, 895–912.
137
Niu, G.-Y, Z.-L. Yang, K. E. Mitchell, F. Chen, M. B. Ek, M. Barlage, A. Kumar, K. Manning,
D. Niyogi, E. Rosero, M. Tewari, Y. Xia, 2011: The community Noah land surface model
with multiparameterization options (Noah-MP): 1. Model description and evaluation with
local-scale measurements. J. Geophys. Res., 116, D12109.
Noh, Y., W.G. Cheon, S.-Y. Hong, and S. Raasch, 2003: Improvement of the K-profile model
for the planetary boundary layer based on large eddy simulation data. Bound.-Layer Meteor.,
107, 401–427.
Noilhan, J., and S. Planton, 1989: A simple parameterization of land surface processes for
meteorological models. Mon. Wea. Rev., 117, 536–549.
Oleson, K. W., et al., 2010: Technical description of version 4 of the Community Land Model
(CLM). NCAR Tech. Note NCAR/TN-478+STR. 266 pp.
Ooyama K. V., 1990: A thermodynamic foundation for modeling the moist atmosphere, J.
Atmos. Sci., 47, 2580–2593.
Pan, H.-L., and W.-S. Wu, 1995: Implementing a mass flux convective parameterization package
for the NMC Medium-Range Forecast model. NMC Office Note 409, 40 pp.
Park, Sungsu, and Christopher S. Bretherton, 2009: The University of Washington shallow
convection and moist turbulence schemes and their impact on climate simulations with the
Community Atmosphere Model. J. Climate, 22, 34493469. doi:10.1175/2008JCLI2557.1
Park, S.-H., W. Skamarock, J. Klemp, L. Fowler, and M. Duda, 2013: Evaluation of global
atmospheric solvers using extensions of the Jablonowski and Williamson baroclinic wave test
case. Mon. Wea. Rev., 141, 3116–3129.
Parrish, D. F., and J. C. Derber, 1992: The National Meteorological Center’s Spectral Statistical
Interpolation analysis system. Mon. Wea. Rev., 120, 1747–1763.
Paulson, C. A., 1970: The mathematical representation of wind speed and temperature profiles
in the unstable atmospheric surface layer. J. Appl. Meteor., 9, 857–861.
Pleim, J. E., 2006: A simple, efficient solution of flux-profile relationships in the atmospheric
surface layer, J. Appl. Meteor. and Clim., 45, 341–347
138
Pleim, J. E., 2007: A combined local and non-local closure model for the atmospheric boundary
layer. Part 1: Model description and testing, J. Appl. Meteor. and Climatol., 46, 1383–1395.
Pleim, J.E. and R.C. Gilliam, 2008: An indirect data assimilation scheme for deep soil temper-
ature in the Pleim-Xiu land surface model. J. Appl. Meteor. Climatol., 48, 1362–1376.
Pleim, J. E. and A. Xiu, 1995: Development and testing of a surface flux and planetary boundary
layer model for application in mesoscale models. J. Appl. Meteor., 34, 16–32.
Pleim, J. E., and A. Xiu, 2003: Development of a land surface model. Part II: Data Assimilation.
J. Appl. Meteor., 42, 1811–1822.
Pollard, R. T., P. B. Rhines and R. O. R. Y. Thompson, 1973: The deepening of the wind-mixed
layer. Geophys. Fluid Dyn., 3, 381–404.
Price, J. F., T. B. Sanford, and G. Z. Forristall, 1994: Forced stage response to a moving
hurricane. J. Phy. Oceanogr., 24, 233–260.
Purser, R. J., W. -S. Wu, D. F. Parrish, and N. M. Roberts, 2003: Numerical aspects of the
application of recursive filters to variational statistical analysis. Part I: Spatially homogeneous
and isotropic Gaussian covariances. Mon. Wea. Rev., 131, 1524–1535.
Rabier, F., H., J. N. Thepaut, and P. Courtier, 1998: Extended assimilation and forecast
experiments with a four-dimensional variational assimilation system. Quart. J. Roy. Meteor.
Soc., 124, 1861–1887.
Rodgers, C. D., 1968: Some extensions and applications of the new random model for molecular
band transmission. Quart. J. Roy. Meteor. Soc., 94, 99–102.
Rutledge, S. A., and P. V. Hobbs, 1984: The mesoscale and microscale structure and organization
of clouds and precipitation in midlatitude cyclones. XII: A diagnostic modeling study of
precipitation development in narrow cloud-frontal rainbands. J. Atmos. Sci., 20, 2949–2972.
Ryan, B. F., 1996: On the global variation of precipitating layer clouds. Bull. Amer. Meteor.
Soc., 77, 53–70.
139
Salamanca, F., and A. Martilli, 2010: A new building energy model coupled with an urban
canopy parameterization for urban climate simulations – part II. Validation with one dimen-
sion off-line simulations. Theor. Appl. Climatol., 99, 345–356.
Sasamori, T., J. London, and D. V. Hoyt, 1972: Radiation budget of the Southern Hemisphere.
Meteor. Monogr., 13, No. 35, 9–23.
Schwartz, C. S., Z. Liu, X.-Y. Huang, 2015: Sensitivity of Limited-Area Hybrid Variational-
Ensemble Analyses and Forecasts to Ensemble Perturbation Resolution. Mon. Wea. Rev.,
143, 3454–3477.
Schwarzkopf, M. D., and S. B. Fels, 1985: Improvements to the algorithm for computing CO2
transmissivities and cooling rates. J. Geophys. Res., 90 (ND6), 541–550.
Schwarzkopf, M. D., and S. B. Fels, 1991: The simplified exchange method revisited — An
accurate, rapid method for computation of infrared cooling rates and fluxes. J. Geophys.
Res., 96 (D5), 9075–9096.
Shin, H. H., and S.-Y. Hong, 2015: Representation of the subgrid-scale turbulent transport in
convective boundary layers at gray-zone resolutions. Mon. Wea. Rev., 143, 250–271.
Shutts, G. J., 2005: A kinetic energy backscatter algorithm for use in ensemble prediction
systems. Quart. J. Roy. Meteor. Soc., 612, 3079–3102.
Skamarock W. C. and J. B. Klemp, 1992: The Stability of Time-Split Numerical Methods for
the Hydrostatic and the Nonhydrostatic Elastic Equations, Mon. Wea. Rev.: 120, 2109–2127.
Smirnova, T. G., J. M. Brown, and S. G. Benjamin, 1997: Performance of different soil model
configurations in simulating ground surface temperature and surface fluxes. Mon. Wea. Rev.,
125, 1870–1884.
140
Stauffer D. R., and N. L. Seaman, 1990: Use of four-dimensional data assimilation in a limited-
area mesoscale model. Part I: Experiments with synoptic-scale data. Mon. Wea. Rev., 118,
1250–1277.
Stauffer D. R., and N. L. Seaman, 1994: Multiscale four-dimensional data assimilation. J. Appl.
Meteor., 33, 416–434.
Stephens, G. L., 1978: Radiation profiles in extended water clouds. Part II: Parameterization
schemes, J. Atmos. Sci., 35, 2123–2132.
Subin, Z.M., Riley, W.J. and Mironov, D. 2012. Improved lake model for climate simulations,
J. Adv. Model. Earth Syst., 4, M02001.
Sukoriansky, S., B. Galperin, and V. Perov, 2005: Application of a new spectral model of
stratified turbulence to the atmospheric boundary layer over sea ice. Bound.-Layer Meteor.,
117, 231–257.
Sun, J., and H. Wang, 2013: Radar Data Assimilation with WRF 4D-Var. Part II: Comparison
with 3D-Var for a Squall Line over the U.S. Great Plains. Mon. Wea. Rev., 141, 2245?-2264.
Sun, J., H. Wang, W. Tong, Y. Zhang, C. Lin, and D. Xu, 2016: Comparison of the Impacts of
Momentum Control Variables on High-Resolution Variational Data Assimilation and Precip-
itation Forecasting. Mon. Wea. Rev., 144, 149–169.
Sun, S., and Y. Xue, 2001: Implementing a new snow scheme in Simplified Simple Biosphere
Model (SSiB), Adv. Atmos. Sci., 18, 335–354.
Tao, W.-K., and J. Simpson, 1993: The Goddard cumulus ensemble model. Part I: Model
description. Terr. Atmos. Oceanic Sci., 4, 35–72.
Tao, W.-K,, J. Simpson, D. Baker, S. Braun, M.-D. Chou, B. Ferrier, D. Johnson, A. Khain, S.
Lang, B. Lynn, C.-L. Shie, D. Starr, C.-H. Sui, Y. Wang, and P. Wetzel, 2003: Microphysics,
radiation and surface processes in the Goddard Cumulus Ensemble (GCE) model. Meteor.
and Atmos. Phys., 82, 97–137.
Tao, W.-K., J. Simpson, and M. McCumber 1989: An ice-water saturation adjustment, Mon.
Wea. Rev., 117, 231–235.
Tao, W.-K., D. Wu, S. Lang, J.-D. Chern, C. Peters-Lidard, A. Fridlind, and T. Matsui (2016),
High-resolution NU-WRF simulations of a deep convective-precipitation system during MC3E:
Further improvements and comparisons between Goddard microphysics schemes and obser-
vations. J. Geophys. Res. Atmos., 121, 1278–1305, doi:10.1002/2015JD023986.
Tiedtke, M., 1989: A comprehensive mass flux scheme for cumulus parameter-
ization in largescale models. Mon. Wea. Rev., 117, 17791800. doi:10.1175/1520-
0493(1989)117¡1779:ACMFSF¿2.0.CO;2.
Thompson, G., R. M. Rasmussen, and K. Manning, 2004: Explicit forecasts of winter precipita-
tion using an improved bulk microphysics scheme. Part I: Description and sensitivity analysis.
Mon. Wea. Rev., 132, 519–542.
141
Thompson, G., P. R. Field, R. M. Rasmussen, and W. D. Hall, 2008: Explicit Forecasts of
Winter Precipitation Using an Improved Bulk Microphysics Scheme. Part II: Implementation
of a New Snow Parameterization. Mon. Wea. Rev., 136, 5095–5115.
Thompson, G., and T. Eidhammer, 2014: A study of aerosol impacts on clouds and precipitation
development in a large winter cyclone. J. Atmos. Sci., 71, 3636–3658.
Tong, W., G. Li, J. Sun, X. Tang and Y. Zhang, 2016: Design Strategies of an Hourly Up-
date 3DVAR Data Assimilation System for Improved Convective Forecasting. Weather and
Forecasting, 31, 1673–1695.
Vendrasco, E. P., J. Sun, D. L. Herdies, and C. F. de Angelis, 2016: Constraining a 3DVAR
Radar Data Assimilation System with Large-Scale Analysis to Improve Short-Range Precipi-
tation Forecasts. J. Appl. Meteor. Climatol., 55, 673–690.
Walko, R. L., W. R. Cotton, M. P. Meyers, and J. Y. Harrington, 1995: Newe RAMS cloud
microphysics parameterization. Part I: The single-moment scheme. Atmos. Res., 38, 29–62.
Wang, W., D. Barker, C. Bruyère, M. Duda, J. Dudhia, D. Gill, J. Micha-
lakes, and S. Rizvi, 2008: WRF Version 3 Modeling System User’s Guide.
http://www.mmm.ucar.edu/wrf/users/docs/user guide V3/.
Wang, H., Skamarock, W. C., Feingold, G., 2009: Evaluation of scalar advection schemes in the
Advanced Research WRF model using large-eddy simulations of aerosol-cloud interactions.
Mon. Wea. Rev., 137, 2547–2558.
Wang, H., J. Sun, S. Fan, and X.-Y. Huang, 2013: Indirect Assimilation of Radar Reflectivity
with WRF 3D-Var and Its Impact on Prediction of Four Summertime Convective Events. J.
Appl. Meteor. Climatol., 52, 889–902.
Wang, H., J. Sun, X. Zhang, X.-Y. Huang, and T. Auligne, 2013: Radar data assimilation with
WRF 4D-Var: Part I. System development and preliminary testing. Mon. Wea. Rev., 141,
2224–2244.
Wang, X., D. M. Barker, C. Snyder, and T. M. Hamill, 2008: A hybrid ETKF-3DVAR data
assimilation scheme for the WRF model. Part I: Observing system simulation experiment.
Mon. Wea. Rev., 136, 5116–5131.
Wang, X., D. M. Barker, C. Snyder, and T. M. Hamill, 2008: A Hybrid ETKF-3DVAR Data
Assimilation Scheme for the WRF Model. Part II: Real Observation Experiments. Mon. Wea.
Rev., 136, 5132–5147.
Webb, E. K., 1970: Profile relationships: The log-linear range, and extension to strong stability,
Quart. J. Roy. Meteor. Soc., 96, 67–90.
Wicker, L. J. and W. C. Skamarock, 2002: Time splitting methods for elastic models using
forward time schemes, Mon. Wea. Rev., 130, 2088–2097.
Wicker, L. J., and R. B. Wilhelmson, 1995: Simulation and analysis of tornado development
and decay within a three-dimensional supercell thunderstorm. J. Atmos. Sci., 52, 2675–2703.
142
Wilson, T. H., and R. G. Fovell, 2018: Modeling the evolution and life cycle of radiative cold
pools and fog. Weather and Forecasting. 33, 203–220.
Wu, W. -S., R. J. Purser, and D. F. Parrish, 2002: Three-Dimensional Variational Analysis with
Spatially Inhomogeneous Covariances. Mon. Wea. Rev., 130, 2905–2916.
Xiao, Q. N., Y. H. Kuo, J. Sun, W. C. Lee, E. Lim, Y. R. Guo, and D. M. Barker, 2005:
Assimilation of Doppler Radar Observations with a Regional 3D–Var System: Impact of
Doppler Velocities on Forecasts of a Heavy Rainfall Case. J. Appl. Met., 44(6), 768–788.
Xiao, Q., Ying-Hwa Kuo, Juanzhen Sun, Wen-Chau Lee, D. M. Barker, and Eunha Lim, 2007:
An approach of radar reflectivity data assimilation and its assessment with the inland QPF
of Typhoon Rusa (2002) at landfall. J. Appl. Meteor. Climat., 46, 14–22.
Xiao, Q., J. Sun, 2007: Multiple radar data assimilation and short-range quantitative precipita-
tion forecasting of a squall line observed during IHOP 2002. Mon. Wea. Rev., 135, 3381–3404.
Xiao, Q., Eunha Lim, D.-J. Won, J. Sun, W.-C. Lee, M.-S. Lee, W.-J. Lee, J.-Y. Cjo, Y.-H.
Kuo, D. M. Barker, D.-K. Lee, and H.-S. Lee, 2008: Doppler radar data assimilation in KMAs
operational forecasting. Bull. Amer. Meteor. Soc., 89, 39–43.
Xiu, A. and J. E. Pleim, 2001: Development of a land surface model part I: Application in a
mesoscale meteorology model. J. Appl. Meteor., 40, 192–209.
Xu, D., Z. Liu, X.-Y. Huang, J. Min, and H. Wang, 2013: Impact of Assimilating IASI Radiance
Observations on Forecasts of Two Tropical Cyclones. Meteorology and Atmospheric Physics,
122, 1–18.
Xu, D., Huang, X.-Y., Liu, Z., Min, J. 2014: Comparisons of Two Cloud Detection Schemes for
Infrared Radiance Observations. Atmospheric and Oceanic Science Letters, 7(4), 358–363.
Xu, D., T. D. Auligne, and X. -Y. Huang, 2015: A validation of the multivariate and mini-
mum residual method for cloud retrieval using radiance from multiple satellites. Advances in
Atmospheric Sciences, 32, 349–362.
Xu, M., Y. Liu, C. Davis and T. Warner, 2002: Sensitivity of nudging parameters on the
performance of a mesoscale FDDA system: A case study. 15th Conference on Numerical
Weather Prediction, 12-16 August, 2002, San Antonio, Texas, 127-130.
Xue, M., 2000: High-order monotonic numerical diffusion and smoothing. Mon. Wea. Rev., 128,
2853–2864.
Xue, Y., P. J. Sellers, J. L. Kinter, and J. Shukla, 1991: A simplified biosphere model for global
climate studies. J. Climate, 4, 345–364.
143
Yang, P., L. Bi, B.A. Baum, K.-N. Liou, G.W. Kattawar, M.I. Mishchenko, and B. Cole, 2013:
Spectrally consistent scattering, absorption, and polarization properties of atmospheric ice
crystals at wavelengths from 0.2 to 100 µm. J. Atmos. Sci., 70, 330–347, doi:10.1175/JAS-D-
12-039.1
Yang, C., Z. Liu, J. Bresch, S.R.H. Rizvi, X.-Y. Huang, and J. Min, 2016: AMSR2 all-sky
radiance assimilation and its impact on the analysis and forecast of Hurricane Sandy with a
limited-area data assimilation system. Tellus A, 68, 30917.
Yang, C., Z. Liu, F. Gao, P. P. Childs, and J. Min, 2017: Impact of assimilating GOES imager
clear-sky radiance with a rapid refresh assimilation system for convection-permitting forecast
over Mexico. J. Geophys. Res. Atmos., 122, 5472–5490.
Zalesak, S. T., 1979: Fully multidimensional flux-corrected transport algorithms for fluids. J.
Comp. Phys., 31, 335362.
Zhang, Chunxi, Yuqing Wang, and Kevin Hamilton, 2011: Improved representation of bound-
ary layer clouds over the southeast pacific in ARWWRF using a modified Tiedtke cumulus
parameterization scheme. Mon. Wea. Rev., 139, 34893513. doi:10.1175/MWR-D-10-05091.1.
Zhang, C., and Y. Wang, 2017: Projected Future Changes of Tropical Cyclone Activity over the
Western North and South Pacific in a 20-km-Mesh Regional Climate Model. J. Climate, 30,
5923-5941. doi:10.1175/JCLI-D-16-0597.1.
Zhang, D.-L., and R.A. Anthes, 1982: A high-resolution model of the planetary boundary layer–
sensitivity tests and comparisons with SESAME–79 data. J. Appl. Meteor., 21, 1594–1609.
Zhang, G. J., and N. A. McFarlane, 1995: Sensitivity of climate simulations to the parame-
terization of cumulus convection in the Canadian Climate Centre general circulation model.
Atmos.Ocean, 33, 407446, doi:10.1080/07055900.1995.9649539
Zhang, X., X.-Y. Huang, and N. Pan., 2013: Development of the Upgraded Tangent Linear and
Adjoint of the Weather Research and Forecasting (WRF) Model. J. Atmos. Oceanic Technol.
30, 1180–1188.
Zhang, X., , X.-Y. Huang, J. Liu, J. Poterjoy, Y. Weng, F. Zhang, and H. Wang, 2014: Devel-
opment of an Efficient Regional Four-Dimensional Variational Data Assimilation System for
WRF. J. Atmos. Oceanic Technol., 31, 2777–2794.
Zhang, X.,Y.-H. Kuo, S.-Y. Chen, X.-Y. Huang, and L.F. Hsiao, 2014: Parallelization Strategies
for the GPS Radio Occultation Data Assimilation with a Nonlocal Operator in the Weather
Research and Forecasting Model. J. Atmos. Oceanic Technol., 31, 2008–2014.
144
Zhang, X., H. Wang, X.-Y. Huang, F. Gao, and N. A. Jacobs, 2015: Using Adjoint-Based
Forecast Sensitivity Method to Evaluate TAMDAR Data Impacts on Regional Forecasts.
Adv. Meteorol., 2015, doi:10.1155/2015/427616.
Zheng, Y., K. Alapaty, J.A. Herwehe, A.D. Del Genio, and D. Niyogi, 2016: Improving high-
resolution weather forecasts using the Weather Research and Forecasting (WRF) Model with
an updated Kain-Fritsch scheme. Mon. Wea. Rev., 144, 833-860.
Zilitinkevich, S. S., 1995: Non-local turbulent transport: pollution dispersion aspects of coher-
ent structure of convective flows, Air Pollution III — Volume I. Air Pollution Theory and
Simulation, Eds. H. Power, N. Moussiopoulos and C.A. Brebbia. Computational Mechanics
Publications, Southampton Boston, 53–60.
145