## Abstract

Urban heat island–induced circulation and convection in three dimensions are investigated theoretically and numerically in the context of the response of a stably stratified uniform flow to specified low-level heating that represents an urban heat island. In a linear, theoretical part of the investigation, an analytic solution for the perturbation vertical velocity in a three-dimensional, time-dependent, hydrostatic, nonrotating, inviscid, Boussinesq airflow system is obtained. The solution reveals a typical internal gravity wave field, including low-level upward motion downwind of the heating center. Precipitation enhancement observed downwind of urban areas may be partly due to this downwind upward motion. The comparison of two- and three-dimensional flow fields indicates that the dispersion of gravity wave energy into an additional dimension results in a faster approach to a quasi-steady state and a weaker quasi-steady flow well above the concentrated heating region in three dimensions.

In a nonlinear, numerical modeling part of the investigation, extensive dry and moist simulations using a nonhydrostatic, compressible model with advanced physical parameterizations [Advanced Regional Prediction System (ARPS)] are performed. While the maximum perturbation vertical velocity in the linear internal gravity wave field exists in the downwind region close to the heating center, the maximum updraft in three-dimensional dry simulations propagates downwind and then becomes quasi stationary. In three-dimensional moist simulations, it is demonstrated that the downwind upward motion induced by an urban heat island can initiate moist convection and result in downwind precipitation. The cloud induced by the downwind upward motion grows rapidly to become deep convective clouds. Heavy rainfalls are localized in a region not far from the heating center by a convective precipitating system that is nearly stationary. The differences in results between two and three dimensions are explained by the presence of (moist) convergence in an additional dimension. The numerical simulation results indicate that the intensity and horizontal structure of the urban heat island affect those of circulation and convection and hence the distribution of surface precipitation.

## 1. Introduction

Numerous observational studies have shown that over and downwind of urban areas the frequency of cloud occurrence and lightning is increased and precipitation is enhanced. Orville et al. (2001) analyzed lightning flash data from the National Lightning Detection Network (NLDN) for the 12-yr period 1989–2000 and found an enhancement of cloud-to-ground lightning over and downwind of Houston, Texas. Shepherd et al. (2002) analyzed rainfall rates derived from the Tropical Rainfall Measuring Mission (TRMM) precipitation radar around major metropolitan areas. The results reveal an average increase of about 28% in monthly rainfall rates within 30–60 km downwind of the metropolis. Using the National Oceanic and Atmospheric Administration (NOAA) satellite images for an 11-yr period, Inoue and Kimura (2004) showed that the low-level cloud frequency is higher over the Tokyo, Japan, metropolitan area than over rural areas on clear summer days. From a climatological examination of summer-season daily precipitation data from 1953 to 2002 for a network of 30 stations in the southeastern United States, Diem and Mote (2005) indicated that precipitation within and downwind of urbanized Atlanta, Georgia, may have been enhanced by urban effects. More observational studies of urban-induced precipitation are summarized in a review paper by Shepherd (2005).

The causes for urban-induced or urban-modified convective phenomena are proposed, which include the urban heat island, increased urban surface roughness, and increased cloud condensation nuclei (CCN) concentration due to urban air pollution. From observations taken from the TRMM satellite, Rosenfeld (2000) found evidence that urban and industrial air pollution can suppress rain and snow by increasing the concentration of CCN, which nucleate many small cloud droplets that coalesce inefficiently into raindrops. On the other hand, Andreae et al. (2004) indicated that the reduced cloud droplet size, and therefore the delay of precipitation onset, can allow invigoration of the updrafts, causing intense ice precipitation, hail, and lightning. Previous studies of the role of increased CCN concentration in suppressing or increasing convective activity are not conclusive (Shepherd 2005), and increased CCN concentration in urban areas does not seem to be a universal cause for the convective phenomena. By analyzing surface meteorological data in the Atlanta area, Bornstein and Lin (2000) documented convective thunderstorms that were initiated by a convergence zone induced by the urban heat island. A two-dimensional modeling study of Baik et al. (2001) shows that the urban heat island can initiate an updraft cell downwind of the heat island, which can initiate moist convection when the thermodynamic conditions are favorable. Through numerical simulations, Rozoff et al. (2003) showed that the urban heat island plays the largest role in initiating deep, moist convection downwind of urban areas. They also showed that increased urban surface roughness induces convergence on the upwind side of urban areas, but this convergence is not strong enough to initiate moist convection. Based upon these studies, the urban heat island appears to be a dominant cause for the observed enhancement of downwind precipitation.

Urban heat island–induced circulation can be theoretically investigated in the context of the response of a stably stratified flow to specified low-level or surface heating that represents an urban heat island. Olfe and Lee (1971) considered urban convection effects by carrying out steady, linearized flow calculations for the velocity and temperature perturbations induced by various surface temperature distributions in a uniform flow. The computed flow field exhibits a downward motion upwind and an upward motion downwind of the heat island center. This is in good qualitative agreement with the observations. Following his earlier work on steady convection over heat islands in the absence of basic-state flow (Kimura 1975), Kimura (1976) examined the effects of uniform basic-state flows on heat island convection by solving a steady-state, linear problem. The convective patterns are classified into three types according to the intensity of the basic-state flow. Lin and Smith (1986) solved a time-dependent, linearized problem in a uniform flow analytically for the airflow over a local heat source and explained a curious negative phase relationship between heating and vertical displacement. Baik (1992) examined the steady-state, linear response of a shear flow to low-level heating and indicated that the magnitude of the perturbation vertical velocity in a shear flow is larger than that in a uniform flow. This is because the basic-state wind shear is a source of the perturbation wave energy. Baik and Chun (1997) dealt with a steady-state, weakly nonlinear problem to investigate the effects of nonlinearity on the airflow past an urban heat island. They showed that when the nondimensional heating depth is large, but still within a valid range of the perturbation expansion, the linear and weakly nonlinear effects constructively work together to produce enhanced upward motion on the downstream side.

These previous theoretical studies are limited to two-dimensional airflow systems to simplify the problem mathematically. However, in an effort to understand urban heat island–induced circulation better, it is necessary to consider a three-dimensional airflow system. Although some theoretical studies (Lin 1986; Lin and Li 1988) considered a three-dimensional airflow system, those studies were directed at understanding the dynamics of deep convection. The first objective of this study is to analytically solve a three-dimensional, transient problem in a uniform flow with low-level heating and examine three-dimensional circulation induced by an urban heat island.

Although linear or weakly nonlinear theoretical studies can provide some basic dynamical insight into urban heat island–induced circulation, these are constrained by the linear or weakly nonlinear assumption. Naturally, to further enhance our understanding of urban heat island–induced circulation and convection, a nonlinear numerical modeling study is needed. This study extends our previous two-dimensional numerical work (Baik et al. 2001) and investigates the effects of nonlinearity on three-dimensional circulation induced by an urban heat island. It will be demonstrated that downwind upward motion induced by an urban heat island can indeed initiate moist convection and result in downwind precipitation in three dimensions and that the intensity and horizontal structure of the urban heat island can affect those of circulation and convection and hence the distribution of surface precipitation. This is the second objective of this study. The theoretical part is given in section 2 and the numerical modeling part is given in section 3. The summary and conclusions are given in section 4.

## 2. Theoretical approach

### a. Governing equations and solutions

Consider a three-dimensional, time-dependent, hydrostatic, nonrotating, inviscid, Boussinesq airflow system. The equations governing perturbations with constant basic-state horizontal wind and buoyancy frequency in the presence of diabatic forcing can be expressed by

Here, *u*, *υ*, and *w* are *x*-, *y*-, and *z*-directional perturbation velocities, respectively; *π* is the perturbation kinematic pressure; *b* is the perturbation buoyancy; *U* and *V* are *x*- and *y*-directional basic-state horizontal wind speeds, respectively; *N* is the buoyancy frequency; *g* is the gravitational acceleration; *c _{p}* is the specific heat of air at constant pressure; and

*T*

_{0}is the reference temperature. The diabatic forcing

*q*, which is regarded as representing an urban heat island, is specified to be bell shaped in the horizontal and to decrease exponentially with height:

where *q*_{0} is the heating amplitude; *a _{x}* and

*a*determine the

_{y}*x*- and

*y*-directional horizontal sizes of the heating, respectively;

*h*is the

*e*-folding heating depth; and

*δ*is the Dirac delta function. Here, the Green function method is used to solve the problem. That is, the solution to the pulse forcing is obtained first and then this Green function is integrated with respect to time to obtain the solution to the steady forcing.

Let us introduce the following nondimensional variables:

where *L* is the horizontal length scale of diabatic forcing, *U*_{0} = (*U*^{2} + *V*^{2})^{1/2} is the magnitude of the basic-state wind speed, and the tilde denotes nondimensional variables. Substituting the above nondimensional variables into (1)–(5) gives (in what follows the tildes are dropped)

where Fr = *U*_{0}/*Nh* is the Froude number and *μ* = *gq*_{0}*L*/(*c _{p}T*

_{0}

*N*

^{2}

*U*

_{0}

*h*) is the nonlinearity factor of thermally induced finite-amplitude waves in three dimensions. The above nondimensional equation set indicates that the linear solution (

*μ*→ 0) depends only on the Froude number for a particular geometry of the heat island and the nonlinear solution depends on the Froude number and the nonlinearity factor. In the case that

*z*and

*h*are nondimensionalized by

*N*/

*U*

_{0},

*μ*becomes equal to the nonlinearity factor in Chun (1991) and Lin and Chun (1991).

For a linear airflow system (*μ* → 0), (8)–(12) can be combined to yield a single equation for the perturbation vertical velocity:

Taking the double Fourier transform in *x* (→*k*) and *y* (→*l*) and the Laplace transform in *t* (→*s*) gives

where

The general solution of (14) is

Two unknown coefficients *A*(*k*, *l*, *s*) and *B*(*k*, *l*, *s*) are determined by imposing a flat-bottom boundary condition [*ŵ*(*k*, *l*, 0, *s*) = 0] and an upper radiation condition [*B*(*k*, *l*, *s*) = 0 for *U* > 0 and *V* = 0]. Then, (15) becomes

The solution in physical space can be obtained by taking the inverse Laplace transform and the inverse double Fourier transform. Subsequent procedures to obtain solution for the perturbation vertical velocity and the solution are given in the appendix.

### b. Results and discussion

Figure 1 shows the time evolution of the perturbation vertical velocity field along *y* = 0 in the steady heating case. The specified parameter values are *U* = 1, *V* = 0, Fr = 1.14, and *a _{x}* =

*a*= 1 (i.e.,

_{y}*U*= 4 m s

^{−1},

*V*= 0 m s

^{−1},

*N*= 0.01 s

^{−1},

*h*= 350 m, and

*a*=

_{x}*a*=

_{y}*L*= 10 km). The

*y*-directional component of the basic-state wind is not considered and specified low-level heating has a circular-shaped structure in the horizontal. In this and following figures, the concentrated heating region is located inside the dotted line.

At *t* = 0.024, an upward motion induced by the low-level heating exists in the vertical near the heating center and a weak compensating downward motion is observed on both sides of the upward motion. At *t* = 0.48, the low-level upward motion intensifies and a slight upstream tilt of the upper part of the perturbation vertical velocity field appears. As time progresses, this pattern becomes apparent and a weak downward motion is present above the low-level upward motion near the heating center. At *t* = 1.44, a moving mode, whose center propagates with a speed equal to the basic-state horizontal wind (Baik et al. 1999), is observed in the region centered at *x* = 1.44. As the moving mode propagates downwind, the region of low-level upward motion is extended farther downwind and the upwind downward motion is enhanced. Precipitation enhancement observed downwind of urban areas may be partly due to this downwind upward motion (Lin and Smith 1986; Baik 1992; Baik and Chun 1997). In the stationary mode that develops near the steady heating region, an upstream phase tilt of upward and downward motion, implying upward wave energy propagation, is observed. At *t* = 2.88, these characteristics of the flow response field become clear and the moving and stationary modes are distinctly separated.

Figure 2 shows the time evolution of the perturbation vertical velocity field at a height of *z* = 2. The specified parameter values are the same as in Fig. 1. At the initial stage, the response to the low-level heating takes the form of a circular upward motion near the heating center, where the heating is concentrated, and a weak compensating downward motion over a wide region surrounding the upward motion. As time progresses, the region of intensifying upward motion extends downwind and the maximum perturbation vertical velocity is observed immediately downwind of the heating center. At *t* = 2.88, a V-shaped region of relatively strong upward motion is located within the area of nearly circular upward motion extending downwind. A moving mode of weak downward motion centered near *x* = 2.88 appears ahead of this V-shaped region. At *t* = 4.32, downwind of the V-shaped upward motion, a similarly shaped weak downward motion appears. Using a group velocity argument, Smith (1980) showed that the wave energy is concentrated near the parabola *y*^{2} = *Nzax*/*U* for the flow aloft past an isolated mountain in a steady state. Lin and Li (1988) extended a group velocity calculation of Smith (1980) to include a nonhydrostatic effect on V-shaped clouds and explained the formation of the repeating, damped oscillations of the disturbance by the nonhydrostatic effect. This feature does not appear in this study because the present system is hydrostatic.

Figure 3 shows perturbation vertical velocity fields along *y* = 0 at *t* = 14.4 with various Froude numbers. A larger Froude number implies a stronger basic-state wind speed for the same buoyancy frequency and heating depth. All of the vertical velocity fields exhibit a low-level upward motion on the downwind side and a relatively weak downward motion on the upwind side, but the magnitude of the perturbation vertical velocity increases as the Froude number decreases. As the Froude number increases, the vertical wavelength becomes large and therefore the tilting slope of upstream phase becomes large for the same horizontal length scale of diabatic forcing. This is apparent from the *e ^{iλz}*

^{/Fr}term in (16), which shows that the vertical wavelength is an increasing function of the Froude number.

To investigate the effects of three dimensionality, the flow response field in three dimensions is compared with that in two dimensions. In two dimensions with a heating structure of *q*(*x*, *z*, *t*) = *q*_{0}[*a*^{2}* _{x}*/(

*x*

^{2}+

*a*

^{2}

*)]*

_{x}*e*

^{−}

*(*

^{z}δ*t*), which is the same as that of Lin and Smith (1986) except for the vertical heating profile, the solution for the perturbation vertical velocity to the steady heating is shown to be

Here,

Note that in the present theoretical study the low-level heating in two dimensions has a different horizontal decay rate from that in three dimensions for a possible Fourier transform {*q*_{2D} ∝ [(*x*/*a _{x}*)

^{2}+ 1]

^{−1},

*q*

_{3D}∝ [(

*x*/

*a*)

_{x}^{2}+ (

*y*/

*a*)

_{y}^{2}+ 1]

^{−3/2}}. However, the comparison of a two-dimensional flow response field with a three-dimensional one with a very elongated heating in the

*y*direction {

*q*

_{3D}→ [(

*x*/

*a*)

_{x}^{2}+ 1]

^{−3/2}} shows that they closely resemble each other except for a small difference in the perturbation magnitude. This implies that any difference resulting from considering different horizontal decay rates of the heating in two and three dimensions would be small in a linear system.

Figure 4 shows perturbation vertical velocity fields at *t* = 14.4 in two and three dimensions (along *y* = 0). The perturbation vertical velocity fields in both two and three dimensions reveal a typical internal gravity wave field in response to thermal forcing. These include an upstream phase tilt of the perturbation field near the steady heating region and a low-level upward (downward) motion downwind (upwind) of the heating center. There are also some differences between the two cases. At an early stage, low-level upward motion near the heating center in three dimensions is much stronger than that in two dimensions. However, after some time, the magnitude of the maximum positive vertical velocity in three dimensions becomes comparable to that in two dimensions. The reason for this faster approach to a quasi-steady state in three dimensions is the dispersion of gravity wave energy into an additional dimension. The weaker stationary gravity wave field well above the concentrated heating region in a quasi-steady state can also be explained by the same reason. The compensating downward motion upwind of the heating center in three dimensions is weaker than that in two dimensions.

## 3. Numerical modeling approach

### a. Numerical model and experimental design

The numerical model used in this study is the Advanced Regional Prediction System (ARPS), which is a three-dimensional, nonhydrostatic, compressible model with advanced physical parameterizations (Xue et al. 2000, 2001). To isolate urban heat island effects on urban-induced circulation and convection, other potential factors such as urban surface roughness are not taken into account. A flat surface is assumed and surface processes, radiation processes, and the rotational effect of the earth are ignored. To parameterize subgrid-scale turbulence, a 1.5-order (turbulent kinetic energy) TKE-based closure scheme is used. Both dry and moist simulations are performed. In moist simulations, the Kessler warm-cloud bulk parameterization scheme is used.

As in the theoretical study, the basic-state horizontal wind is taken to be constant with height and the *y*-directional component is not considered. Also, a constant stratification of *N* = 0.01 s^{−1} is considered, with which the basic-state potential temperature profile is determined (the basic-state potential temperature at the surface is specified as 288 K). The basic-state relative humidity is constant from the surface to *z* = 1 km (RH = 85%), decreases linearly with height up to *z* = 11 km, and is constant above that height (RH = 10%). This thermodynamic sounding has a convective available potential energy (CAPE) of 355 J kg^{−1} when the moisture effect is included. The lifting condensation level (LCL) is 966 hPa and the level of free convection (LFC) is 902 hPa. The domain size is 150 km × 100 km × 15 km. The horizontal and vertical grid sizes are 1 km and 150 m, respectively. To minimize wave reflection at boundaries, a lateral radiation boundary condition is used and a damping layer is put in from *z* = 12 km to the model top height. The model is integrated up to *t* = 6 h with a large time step of 2 s and a small time step, which is related to acoustic waves, of 0.5 s. The structure of specified low-level heating centered at *x* = *y* = 50 km is the same as that in the theoretical study.

### b. Results and discussion

Extensive dry and moist simulations with various heating amplitudes, heating structures, and basic-state wind speeds are performed to examine urban heat island–induced circulation and convection in three dimensions. Unless otherwise stated, *a _{x}* =

*a*= 10 km (

_{y}*γ*=

*a*/

_{y}*a*= 1).

_{x}#### 1) Dry simulations

As mentioned previously, the nonlinearity factor of thermally induced finite-amplitude waves in three dimensions can be expressed by *μ* = *gq*_{0}*L*/(*c _{p}T*

_{0}

*N*

^{2}

*U*

_{0}

*h*). Accordingly, the degree of nonlinearity of the flow increases as the heating amplitude increases. For the larger heating amplitude, the intensity of the internal gravity wave field in a linear airflow system increases in proportion to the heating amplitude, but its pattern remains unchanged because the solution for the perturbation vertical velocity is directly proportional to the heating amplitude. To examine the effects of nonlinearity on three-dimensional circulation induced by an urban heat island, numerical simulations with various heating amplitudes are performed. Figure 5 shows vertical velocity fields along

*y*= 50 km (the heating center) at

*t*= 3 h with

*q*

_{0}= 0.2, 0.4, 0.6, 0.8, and 1.0 J kg

^{−1}s

^{−1}and

*U*= 4 m s

^{−1}(

*μ*= 0.51, 1.02, 1.53, 2.04, and 2.55 and Fr = 1.14). In calculating

*μ*, the horizontal size of the heating in the basic-state wind direction

*a*is regarded as the horizontal length scale of diabatic forcing

_{x}*L*. Simulated flow fields at

*t*= 3 h are in a quasi-steady state, in which the magnitude of the perturbation decreases slowly, but its pattern remains similar. Note that the contour intervals are set to be in proportion to the heating amplitude. When the degree of nonlinearity of the flow is small (Fig. 5a), the flow response field overall resembles the linear internal gravity wave field (Fig. 3d). Some differences between the two fields are largely related to dynamic frame differences (e.g., linear, inviscid, hydrostatic versus nonlinear, viscid, nonhydrostatic). However, an intensified maximum vertical velocity is observed farther downwind of the heating center as the heating amplitude increases (Figs. 5a–c). For heating amplitudes larger than 0.6 J kg

^{−1}s

^{−1}, any further increase in the heating amplitude does not produce a farther downwind shift of the location of the maximum upward motion. A weakly organized updraft develops within the region of downwind upward motion as the heating amplitude increases (Figs. 5c–e).

Figure 6a shows the vertical velocity field along *y* = 50 km at *t* = 3 h with *q*_{0} = 0.8 J kg^{−1} s^{−1} and *U* = 4 m s^{−1} (*μ* = 2.04 and Fr = 1.14). The vertical velocity field is generally similar to the linear internal gravity wave field in response to the heating (Fig. 1), including a relatively strong downwind upward motion with an upstream phase tilt and a weak upwind downward motion. However, the maximum vertical velocity in the linear internal gravity wave field exists in the downwind region close to the heating center (Fig. 1), while the maximum updraft in the three-dimensional dry simulation propagates downwind and becomes quasi stationary 17 km downwind of the heating center where strong convergence and warm advection maximum occur (Figs. 6a–d). After some time, the overall intensity of the vertical velocity field weakens because of the dissipation effect (*w*_{max} = 0.41 m s^{−1} at *t* = 2 h and 0.35 m s^{−1} at *t* = 6 h), but its pattern remains similar. The vertical velocity field at a height of *z* = 1 km shows that there exists a V-shaped region of upward motion, similar to that in the theoretical study (Figs. 2 and 6e).

The nonlinearity factor is not only proportional to the heating amplitude, but also inversely proportional to the basic-state wind speed. Therefore, the degree of nonlinearity of the flow increases with decreasing basic-state wind speed for the same heating amplitude. However, note that the Froude number also depends on the basic-state wind speed. Accordingly, as the basic-state wind speed decreases, the nonlinearity factor increases, but on the other hand the Froude number decreases. Figure 7 shows vertical velocity fields along *y* = 50 km at *t* = 3 h with *q*_{0} = 0.8 J kg^{−1} s^{−1} and *U* = 2, 4, 6, and 8 m s^{−1} (*μ* = 4.09, 2.04, 1.36, and 1.02 and Fr = 0.57, 1.14, 1.71, and 2.29). In the case of *U* = 2 m s^{−1} (Fig. 7a), an organized updraft is present within a narrow region of downwind upward motion, which is not far from the heating center. As the basic-state wind speed increases, the maximum upward motion is located farther away from the heating center and has weaker intensity and the region of low-level upward motion spreads out downwind due to the larger horizontal advection effect. In the cases of *U* = 2, 4, 6, and 8 m s^{−1}, the maximum vertical velocity is located 8, 17, 24, and 27 km downwind of the heating center, respectively. The magnitude of the maximum upward motion with *q*_{0} = 0.8 J kg^{−1} s^{−1} and *U* = 8 m s^{−1} (Fig. 7d) is similar to that with *q*_{0} = 0.4 J kg^{−1} s^{−1} and *U* = 4 m s^{−1} (Fig. 5b), since both cases have the same nonlinearity (*μ* = 1.02). However, in the case of *q*_{0} = 0.8 J kg^{−1} s^{−1} and *U* = 8 m s^{−1}, larger advection effect results in a farther downwind movement of the region of the maximum vertical velocity compared with the case of *q*_{0} = 0.4 J kg^{−1} s^{−1} and *U* = 4 m s^{−1}. Extensive dry simulations with various basic-state wind speeds indicated that the location of the maximum updraft in three dimensions becomes quasi stationary after some time regardless of the basic-state wind speed.

The above results in three dimensions are quite different from those in two dimensions showing that a well-organized downwind updraft cell, which is separated from the stationary gravity wave field near the heating center, propagates downwind continuously (Baik et al. 2001). Differences in the flow response field between two- and three-dimensional dry simulations are basically due to the difference in the *y*-directional structure of specified low-level heating. Figure 8 shows vertical velocity fields along *y* = 50 km at *t* = 6 h in dry simulations with *q*_{0} = 0.8 J kg^{−1} s^{−1}, *U* = 4 m s^{−1}, *a _{x}* = 10 km (

*μ*= 2.04 and Fr = 1.14), and different values of

*a*(

_{y}*a*=

_{y}*a*, 2

_{x}*a*, 3

_{x}*a*, and 4

_{x}*a*i.e.,

_{x}*γ*= 1, 2, 3, and 4). As the heating becomes more elongated in the

*y*direction, the maximum upward motion is located farther downwind and the strong updraft is more organized. As the three-dimensional heat source is increasingly elongated in the

*y*direction, that is, the

*y*-directional structure of the low-level heating in three dimensions approaches that in two dimensions, it was found that a downwind updraft cell, separated from the stationary gravity wave field, appears and continues to move downwind.

The dependence of flow on the horizontal structure of the urban heat island is investigated by performing two numerical simulations in which elliptical-shaped heating structures are considered. Figure 9 shows vertical velocity fields along *y* = 50 km at *t* = 3 h in dry simulations with (*a _{x}*,

*a*) = (20 km, 5 km) and (5 km, 20 km). Here

_{y}*q*

_{0}and

*U*are specified as 0.8 J kg

^{−1}s

^{−1}and 4 m s

^{−1}, respectively (

*μ*= 4.09 and 1.02, Fr = 1.14, and

*γ*= 0.25 and 4). Both the circular-shaped heating case of

*a*=

_{x}*a*= 10 km and the elliptical-shaped heating case of

_{y}*a*= 20 km and

_{x}*a*= 5 km or

_{y}*a*= 5 km and

_{x}*a*= 20 km are set to have the same total amount of heating. Note that the total amount of heating considered in this study is proportional to the product of

_{y}*a*and

_{x}*a*. In the case of

_{y}*a*= 4

_{x}*a*(=20 km), there exists a broad region of upward motion not only near and downwind of the heating, but also on the upwind side due to the elongated heating structure in the basic-state wind direction. The maximum updraft moves downwind and then becomes quasi stationary, as in the case of circular-shaped heating, but it is closer to the heating center and has a stronger intensity. The response to the heating elongated perpendicularly to the basic-state wind direction (

_{y}*a*= 4

_{y}*a*= 20 km) shows a narrower region of downwind upward motion and a stronger compensating upwind downward motion compared with the case of

_{x}*a*=

_{x}*a*(=10 km) (Fig. 6a). The linear internal gravity wave–like flow response field is observed near the heating region. In particular, no downwind movement of the maximum vertical velocity occurs and it stays right downwind of the heating center as in the theoretical study, even though the flow is nonlinear (

_{y}*μ*= 1.02). The maximum upward motion, located at a higher altitude, is much weaker than that in the circular-shaped heating case.

The dry simulation results clearly indicate that an organized updraft of relatively strong intensity develops within the region of downwind upward motion as the flow becomes nonlinear. An urban heat island–induced downwind upward motion that strengthens as the heating amplitude increases or the basic-state wind speed decreases (i.e., as the degree of nonlinearity of the flow increases) can initiate moist convection and result in downwind precipitation under favorable thermodynamic conditions. It is expected that any change in the horizontal location of the maximum upward motion, resulting from different urban heat island intensities and structures and basic-state wind speeds, can change the location at which moist convection occurs and the distribution of surface precipitation. This will be examined in section 3b(2).

The turbulent diffusion of heat and momentum may not be neglected in the urban heat island problem. In the real atmosphere, heat sources for urban heat islands are located at and near the surface and heat is transferred in the vertical by advection and turbulent diffusion processes. In the present theoretical approach, the turbulent heat diffusion term (also the turbulent momentum diffusion term) is neglected for mathematical tractability, but the heating has a vertical structure that decays with height. Therefore, we might think that the distributed heat in the vertical (a vertical heating structure) results partly or largely from turbulent heat diffusion process. That is, the turbulent heat diffusion process is implicitly reflected in the specified heating to some extent, even if it is neglected in the thermodynamic energy equation. In the numerical simulations, for consistency the heating is specified as having the same structure as that employed in the theoretical approach. An appropriate method might be to add heating at the lowest model level or surface, making heat be distributed in the vertical by model physics such as turbulence process. This method deserves an investigation.

#### 2) Moist simulatons

Figure 10 shows cloud water mixing ratio, rainwater mixing ratio, vertical velocity, and temperature fields along *y* = 50 km at *t* = 3 h in a moist simulation with *q*_{0} = 1 J kg^{−1} s^{−1} and *U* = 4 m s^{−1} (*μ* = 2.55 and Fr = 1.14). The first cloud water is produced 11 km downwind of the heating center at *t* = 62 min. At this time, the maximum downwind updraft induced by the heating is found at the same location with an intensity of 0.49 m s^{−1}. The coincidence of the horizontal location of the first cloud water formation with that of the maximum downwind updraft implies that the downwind upward motion induced by an urban heat island dynamically initiated moist convection. The first rainwater is produced 15 km downwind of the heating center 25 min after moist convection first occurs. The low-level cloud induced by the downwind upward motion grows rapidly to become deep convective clouds with a slight downstream tilt and heavy rainfalls are localized in a region not far from the heating center by a single-storm-like convective system that develops in a weak or no shear environment (Figs. 10a,b). At the upper level downwind of the convective clouds, stratiform clouds composed of cloud water appear. Before the first cloud water is produced, urban heat island–induced vertical velocity field is similar to that of the corresponding dry simulation (Fig. 5e). However, after the downwind upward motion that is intensified up to a certain intensity initiates moist convection, the resultant vertical velocity field shows a strong updraft within the convective system and a downdraft over a surface precipitation region (Fig. 10c). The temperature within the convective system is higher than that in the surroundings due to the latent heat released by condensing water vapor, while downwind of the convective system the temperature near the surface is lower than that in the surroundings due to the evaporative cooling of falling raindrops (Fig. 10d).

The horizontal distribution of 3-h accumulated surface precipitation amount from *t* = 3 to 6 h in the case corresponding to Fig 10 is shown in Fig. 12a. There is surface precipitation over a wide region from *x* = 59 to 87 km (note that the contour interval is 5 mm), but nearly all of it occurs in a localized region, which has an elongated shape in the basic-state wind direction, by a convective precipitating system that is nearly stationary on the downwind side. The maximum 3-h accumulated precipitation amount is 71 mm and is observed 21 km downwind of the heating center. Very weak precipitation occurs downwind of the primary precipitation zone by stratiform clouds located downwind of the deep convective clouds. The precipitation pattern in three dimensions is different from that in two dimensions. In two dimensions, the amount of surface precipitation is greatly reduced and its area, located farther away from the heating center, exhibits a more widespread distribution (Baik et al. 2001).

To explain differences in the amount and horizontal location of surface precipitation between two- and three-dimensional moist simulations, numerical simulations with various *y*-directional sizes of heating are performed. Figures 11a,b show the time and horizontal location of the first cloud water and rainwater formation as a function of *a _{y}* in moist simulations with

*q*

_{0}= 1 J kg

^{−1}s

^{−1}and

*U*= 4 m s

^{−1}(

*μ*= 2.55, Fr = 1.14, and

*γ*= 1, 2, 3, and 4). As the

*y*-directional decay rate of the heating decreases, the time required for the first cloud water and rainwater formation increases and the horizontal location is placed farther downwind due to reduced (moist) convergence in the

*y*direction. Additional numerical experiments with larger values of

*a*showed that the time and horizontal location of the first cloud water and rainwater formation in three dimensions approach those in two dimensions as the

_{y}*y*-directional size of heating becomes very large. Figure 11c shows the horizontal location of the maximum updraft as a function of

*a*in dry simulations at the time the first cloud water is produced in the corresponding moist simulations. The horizontal location of the maximum downwind updraft in dry simulations is exactly the same as that of the first cloud water formation in the corresponding moist simulations. This is also true in the numerical simulations with different heating amplitudes (Fig. 13).

_{y}Figure 12 shows the horizontal distribution of 3-h accumulated surface precipitation amount from *t* = 3 to 6 h with various *y*-directional sizes of heating. As stated above, surface precipitation concentrates in a narrow area not far from the heating region in the case of *a _{y}* =

*a*= 10 km. However, the precipitating region spreads over a broad area farther downwind with increasing

_{x}*a*. In the case of

_{y}*a*= 2

_{y}*a*= 20 km, surface precipitation is observed in a narrow region elongated in the basic-state wind direction as in the circular-shaped heating case, but it exists farther downwind of the heating center compared with that of the circular-shaped heating case (Fig. 12b). In the case of

_{x}*a*= 4

_{y}*a*= 40 km, the precipitating region not only propagates in the downwind direction, but it also occurs over a widespread area in the

_{x}*y*direction (Fig. 12d). Secondary maxima in precipitation amount are observed at locations of (

*x*,

*y*) = (88 km, 41 km) and (87 km, 58 km).

The horizontal sizes of heating in the experiments with *a _{x}* = 20 km,

*a*= 20 km,

_{y}*a*= 30 km, and

_{y}*a*= 40 km are large compared with those of present-day urban heat islands. Those experiments are, however, considered to systematically examine the effects of the horizontal structure of the urban heat island on its induced circulation and convection. The larger horizontal sizes of heating might be thought of as those pertaining to a future situation in which cities have grown to cover a much larger area. The experiments with large heating amplitudes (e.g., Fig. 13) are similarly considered (a systematic examination of the effects of the heating amplitude and an assumption of a future situation with enhanced urban heat islands).

_{y}Our previous two-dimensional moist simulations (Baik et al. 2001) show that as the heating amplitude increases, the first cloud water and rainwater are produced earlier and closer to the heating center. The same feature appears in the present three-dimensional study (Fig. 13), but as mentioned above, the time required for the first cloud water and rainwater formation shortens and the horizontal location is placed closer to the heating center compared with the two-dimensional moist simulations for a given heating amplitude. This is because in three dimensions *y*-directional convergence makes the downwind updraft develop more rapidly at an early stage. The heavy precipitation region is also located closer to the heating center as the heating amplitude increases (not shown). In the cases of *q*_{0} = 0.8, 1.0, 1.2, 1.4, and 1.6 J kg^{−1} s^{−1}, the maximum accumulated precipitation amounts are observed 30, 21, 16, 13, and 10 km downwind of the heating center, respectively. The decreasing basic-state wind speed has the same effect on the cloud water and rainwater formation as increasing heating amplitude. That is, as the basic-state wind speed decreases, the first cloud water and rainwater are produced earlier and closer to the heating center and hence the precipitation region is located closer to the heating center (not shown).

To investigate the influence of the horizontal structure of the urban heat island on moist convection and its resulting precipitation, moist simulations with two different types of elliptical-shaped heating are performed. The *y*-directional width of precipitation region for the elliptical-shaped heating case of *a _{x}* = 20 km and

*a*= 5 km is similar to that in the case of

_{y}*a*=

_{x}*a*= 10 km, but the precipitation region has a more elongated shape in the basic-state wind direction (Figs. 14a,b). As the heating amplitude increases, there is an upwind movement of the precipitation region (not shown), resulting in a large amount of precipitation not only downwind, but also over the urban heat island, unlike the circular-shaped heating cases in which heavy precipitation occurs only on the downwind side. When the heat source is elongated perpendicular to the basic-state wind direction, the precipitation region is observed farther downwind of the heating center compared with the case of circular-shaped heating (Figs. 14a,c).

_{y}The main conclusions in the moist simulation part are that the downwind updraft induced by an urban heat island can initiate moist convection in three dimensions and that the intensity and horizontal structure of the urban heat island affect those of circulation and convection and hence the surface precipitation distribution. The types or structures of simulated convective systems depend on many factors, including specified heating structures and basic-state wind and thermodynamic profiles. In the real atmosphere, the heating would more or less always be present and the basic-state environment varies with time. In this study, we have taken a rather simple situation to get the basic understandings of urban heat island–induced circulation and convection with steady heating and constant basic-state wind and buoyancy frequency. Note again that the heating considered in the present numerical simulations is not a pulse heating but a steady heating as in the theoretical study, but the heating is suddenly imposed at *t* = 0. To examine whether the sudden introduction of the heating at *t* = 0 changes the main conclusions, numerical simulations were performed in which the heating is specified such that *q* is 0 at *t* = 0, gradually increases, and then reaches a steady state (i.e., the heating is transient). It was found that the main conclusions in the moist simulation part are still valid. As should be expected, moist convection starts at a later time (initiation time depends on the transient heating profile).

## 4. Summary and conclusions

In this study, three-dimensional circulation and convection induced by an urban heat island were investigated using a linear theory and a nonlinear numerical model. Our theoretical study was accomplished by examining the three-dimensional, transient response of a stably stratified atmosphere to prescribed low-level heating. The airflow system considered is a three-dimensional, time-dependent, hydrostatic, nonrotating, inviscid, Boussinesq system. The thermal forcing that represents an urban heat island is bell shaped in the horizontal and decreases exponentially with height. The equations governing small-amplitude perturbations in constant basic-state horizontal wind and buoyancy frequency with low-level heating were analytically solved using the Fourier transform method in horizontal space and the Green function method in time. The analytic solution for the perturbation vertical velocity reveals a typical flow response field of thermally induced internal gravity waves, in which the low-level upward motion induced by the heating intensifies and extends downwind with time. This low-level upward motion on the downwind side explains to some extent the observed precipitation enhancement downwind of urban areas. The comparison of two- and three-dimensional flow fields indicates that the dispersion of gravity wave energy into an additional dimension results in a faster approach to a quasi-steady state and a weaker quasi-steady flow above the concentrated heating region in three dimensions.

In a numerical approach using a three-dimensional model with advanced physical parameterizations (ARPS), extensive dry and moist simulations with various urban heat island intensities and horizontal structures and basic-state wind speeds were performed. The nonlinear flow response field in dry simulations shows no significant difference compared with the linear internal gravity wave field except for the horizontal location of the maximum upward motion. Unlike the linear system, in which the maximum upward motion exists in the downwind region close to the heating center, the maximum updraft in the three-dimensional nonlinear system moves in the downwind direction and then becomes quasi stationary. In moist simulations, urban heat island–induced downwind upward motion that is intensified up to a certain intensity initiates moist convection and results in downwind precipitation. Heavy precipitation occurs in a localized region, which has an elongated shape in the basic-state wind direction not far from the heating center, by a convective precipitating system that is nearly stationary on the downwind side. From the numerical simulations with various *y*-directional sizes of heating, different flow responses to the low-level heating between two and three dimensions were explained by the difference in the *y*-directional structure of specified low-level heating. The numerical simulation results indicate that the intensity and horizontal structure of the urban heat island affect those of circulation and convection and hence the distribution of surface precipitation.

The urban heat island is a common phenomenon observed in urban areas. Our theoretical and numerical study highlights a universal role that the urban heat island plays in downwind precipitation enhancement. The current study will be extended to examine the extent to which increased surface roughness contributes to urban-induced circulation and convection, especially to downwind precipitation enhancement, although its role is not thought to dominate over the role that the urban heat island plays. Also, the roles of increased CCN concentration will be examined using a cloud model. Urbanization will be continued and, accordingly, urban-induced or urban-modified weather and climate would be more prominent in the future than at present. More attention to this problem is required.

## Acknowledgments

The authors are very grateful to two anonymous reviewers for providing valuable comments on this work. This work was funded by the Korea Meteorological Administration Research and Development Program under Grant CATER 2006-2202.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

### APPENDIX

#### Analytic Solution for the Perturbation Vertical Velocity

To obtain the solution in physical space, the inverse Laplace transform in *s* (→*t*) and the inverse double Fourier transform in *k* (→*x*) and *l* (→*y*) of (16) are taken. After some manipulation, one can get

where

The solution for the perturbation vertical velocity to the pulse forcing is finally expressed by

Here,

The solution for the perturbation vertical velocity to the steady forcing is obtained by integrating (A2) with respect to time:

## Footnotes

*Corresponding author address*: Jong-Jin Baik, School of Earth and Environmental Sciences, Seoul National University, Seoul 151-742, South Korea. Email: jjbaik@snu.ac.kr