## Abstract

This paper presents high-resolution numerical solutions of a nonlinear zonally symmetric slab model of the intertropical convergence zone (ITCZ) boundary layer. The boundary layer zonal and meridional flows are forced by a specified pressure field, which can also be interpreted as a specified geostrophically balanced zonal wind field *u*_{g}(*y*). One narrow on-equatorial peak in boundary layer pumping is produced when the forcing is easterly geostrophic flow along the equator and two narrow peaks in boundary layer pumping are produced on opposite sides of the equator (a double ITCZ) when the forcing is westerly geostrophic flow along the equator. In the case when easterlies are surrounding a westerly wind burst, once again a double ITCZ is produced, but the ITCZs have significantly more intense boundary layer pumping than the case of only westerly geostrophic flow. A comparison of the numerical solutions to those of classical Ekman theory suggests that the meridional advection term *υ*(∂*υ*/∂*y*) plays a vital role in strengthening and narrowing boundary layer pumping regions while weakening and broadening boundary layer suction regions.

## 1. Introduction

Figure 1 is a GOES visible–infrared blended image of the Pacific Ocean on 11 March 2015, a day when there were twin tropical cyclones in the west and a double intertropical convergence zone (ITCZ) in the east. A striking feature of this image, and many other similar images, is the narrowness of the double ITCZ bands and, hence, the narrowness of the rising parts of the Hadley circulation. The purpose of this paper is to better understand the boundary layer dynamics associated with these ITCZ features. The meridional distribution of Ekman pumping in and around the ITCZ will be examined using a zonally symmetric slab boundary layer model that includes the meridional advection term *υ*(∂*υ*/∂*y*), where *υ* is the meridional velocity. With the inclusion of the *υ*(∂*υ*/∂*y*) term in the boundary layer dynamics, very sharp meridional gradients can appear in the *υ* field, and the Ekman pumping can become very localized. An analogous process in the hurricane eyewall has recently been studied by Smith and Montgomery (2008), Williams et al. (2013), Slocum et al. (2014), and Williams (2015).

For a review of early research on the formation of the ITCZ, the reader is referred to Charney (1969, 1971), Yamasaki (1971), Pike (1971, 1972), Mahrt (1972a,b), Holton et al. (1971), Chang (1973), and Holton (1975). Charney (1969, 1971) suggested that the mechanism that leads to the formation of tropical cyclones might also be responsible for the formation of the ITCZ itself. To support this hypothesis, he presented a linear instability argument for the existence of the ITCZ, with essentially the same physical mechanisms as the Charney and Eliassen (1964) argument for the existence of tropical cyclones, but with line symmetry replacing axisymmetry. Holton et al. (1971) and Holton (1975) approached the topic by using a boundary layer model on the equatorial *β* plane forced by a specified pressure field. The former linearized the boundary layer equations; the latter added the effects of horizontal advection. Both studies found that convergence is frictionally driven and concentrated at critical latitudes, where the disturbance frequency is equal to the Coriolis parameter. Holton (1975) discussed the importance of horizontal advection in strengthening and concentrating regions of convergence while weakening and spreading divergence regions. Despite a general agreement between many of these studies that boundary layer convergence is often concentrated at critical latitudes, these ideas appear to have shifted in the early to mid-1980s, likely because of improved understanding of the role of thermodynamical processes at low levels.

For example, Lindzen and Nigam (1987), Waliser and Somerville (1994), Stevens et al. (2002), McGauley et al. (2004), Sobel and Neelin (2006), Raymond et al. (2006), and Back and Bretherton (2009) illustrated the importance of sea surface temperature (SST) gradients in supporting large-scale pressure gradients, which enhance low-level horizontal wind convergence. The extended Ekman model that does not include horizontal advection (instead, it is a local balance between the Coriolis, pressure gradient, and frictional forces) was used by Stevens et al. (2002), McGauley et al. (2004), Raymond et al. (2006), and Back and Bretherton (2009). Raymond et al. (2006) argued that a scale analysis shows that the *υ*(∂*u*/∂*y*) term is not important. However, it is important to note that a single horizontal scale can often be difficult to determine when multiple-scale features are present such as those observed in the ITCZ boundary layer. A more general model that includes horizontal advection was used by Tomas et al. (1999), who suggest that the *υ*(∂*u*/∂*y*) term is vital to correctly simulate the ITCZ. Sobel and Neelin (2006) also studied the ITCZ boundary layer and free troposphere using an equatorial *β*-plane model that includes both horizontal advection and moisture, as they were specifically interested in narrow ITCZs. They analyze the magnitude of the individual model terms noting that SST gradients associated with the Lindzen and Nigam (1987) effect are most important in determining ITCZ width and intensity. Also, horizontal advection plays an important role near narrow boundary layer convergence regions, as seen in their Fig. 5.

In this paper, we explore the roles of horizontal advection and horizontal diffusion in determining the location of the Ekman pumping. As global models increase their horizontal resolution, nonlinear dynamical effects, such as horizontal advection, may become more important. Also, such nonlinear dynamical effects may have significant consequences for the thermodynamical fields. Therefore, although we realize that the ITCZ is a complex phenomenon, we have chosen to analyze it in an idealized framework focusing on its dynamical aspects. We present a time-dependent dynamical view of the ITCZ, where the pressure gradient in and just above the boundary layer supports both the evolution of the boundary layer zonal and meridional winds along with the formation and narrowing of Ekman pumping and vorticity in the ITCZ due to the nonlinear effects within the meridional momentum equation. The boundary layer we consider is the subcloud boundary layer, where the dynamical and thermodynamical variables tend to be well mixed (Johnson et al. 2001; Bond 1992; Yin and Albrecht 2000; McGauley et al. 2004). Therefore, we do not explicitly treat the influence of SSTs on the boundary layer pressure gradient as done in many of the studies discussed above. We also provide new insight into the topic of the ITCZ boundary layer through our mathematical analysis of the boundary layer equations as a nonlinear hyperbolic system.

This paper is organized in the following manner. Section 2 introduces the slab boundary layer model. Section 3 contains a proof that, in the absence of the horizontal diffusion terms, the slab boundary layer equations constitute a nonlinear hyperbolic system that can therefore be written in characteristic form. This characteristic form is the result of converting the original partial differential equations into ordinary differential equations along characteristics. The classification of the slab model as a nonlinear hyperbolic system alerts us to the possibility of shocks or shocklike structures (Burgers 1948). We define a shocklike structure as a region of rapid changes in the fluid velocity that would be a true discontinuity, or shock, in the absence of horizontal diffusion. Section 4 presents the slab version of classical Ekman theory (Ekman 1905) that can be regarded as the local (i.e., horizontal advection is neglected), steady-state version of the complete nonlinear slab model. In section 5, we perform three experiments using both the fully nonlinear numerical model and classical Ekman theory. We are then able to better assess the important role of horizontal advection in ITCZ formation. Some concluding remarks are presented in section 6, including the implications of the present work on understanding single and double ITCZs.

## 2. Slab boundary layer model

The model considers zonally symmetric, boundary layer motions of an incompressible fluid on the equatorial *β* plane. The frictional boundary layer is assumed to have constant depth *h*, with zonal and meridional velocities *u*(*y*, *t*) and *υ*(*y*, *t*) that are independent of height between the top of a thin surface layer and height *h*, and with vertical velocity *w*(*y*, *t*) at height *h*. In the overlying layer, the meridional velocity is assumed to be negligible and the zonal velocity *u*_{g}(*y*) is assumed to be in geostrophic balance and to be a specified function of *y*. The boundary layer flow is driven by the same meridional pressure gradient force that occurs in the overlying fluid, so that in the meridional equation of boundary layer motion, the pressure gradient force can be expressed as the specified function *βyu*_{g}. The governing system of differential equations for the boundary layer variables *u*(*y*, *t*), *υ*(*y*, *t*), and *w*(*y*, *t*) then takes the form

where

is the wind speed at 10-m height (Powell et al. 2003), *β* = 2Ω/*a*, Ω and *a* are Earth’s rotation rate and radius, and *K* is the constant horizontal diffusivity. The drag factor *c*_{D}*U* is assumed to depend on the 10-m wind speed according to

where the 10-m wind speed *U* is expressed in m s^{−1}, and where (5) applies for *U* ≤ 25 m s^{−1} (Large et al. 1994). The boundary conditions are

The initial conditions are

See the appendix for a complete derivation of (1)–(3) from conservation principles such as absolute angular momentum.

In the absence of the horizontal diffusion terms, the slab boundary layer equations constitute a hyperbolic system that can be written in characteristic form, as presented in the next section. Knowledge of the characteristic form is useful in understanding the possible formation of shocks and shocklike structures.

## 3. Characteristic form

Equations (1)–(7) constitute a nonlinear system with the nonlinearity arising from the *υ*(∂*u*/∂*y*) and *υ*(∂*υ*/∂*y*) terms, the *w* terms, and the *c*_{D}*U* terms. The *υ*(∂*u*/∂*y*) and *υ*(∂*υ*/∂*y*) terms are often referred to as “quasi linear” because they are linear in the first derivatives but the coefficient of these derivatives involves the dependent variable *υ*. Neglecting the horizontal diffusion terms, the system of equations constitutes a nonlinear hyperbolic system; that is, it can be rewritten in characteristic form. An important aspect of nonlinear hyperbolic equations is the possibility of shock formation. Knowledge of the characteristic form allows for a deeper understanding of the way that characteristics can intersect and thereby produce a discontinuity in *υ* and a singularity in *w*. The characteristic form can be derived by rearranging (1) and (2) in such a way that all of the terms involving the derivatives (∂*u*/∂*t*), (∂*u*/∂*y*), (∂*υ*/∂*t*), and (∂*υ*/∂*y*) appear on the left-hand sides and the other terms appear on the right-hand sides. In regions where *w* ≥ 0, the *w* terms in (1) and (2) vanish. In regions where *w* < 0, the *w* terms do not vanish, in which case these terms need to be expressed in terms of (∂*υ*/∂*y*). This procedure is easily accomplished by noting that the mass continuity equation in (3) allows (1) and (2) to be written in the form

where

The forms (8) and (9) are convenient because the nonlinearities associated with spatial derivatives are on the left-hand side, while all of the other linear and nonlinear terms are on the right-hand side. As discussed by Whitham (1974), the classification of the system (8) and (9) as a hyperbolic system and the determination of the characteristic form of this system depends on finding the eigenvalues and left eigenvectors of the matrix

This matrix is defined by the coefficients of the (∂*u*/∂*y*) and (∂*υ*/∂*y*) terms on the left-hand sides of (8) and (9). For *n* = 1, 2, let be the left eigenvector of corresponding to the eigenvalue *λ*^{(n)}; that is,

As can be shown through direct substitution into (13), the two eigenvalues and the two corresponding left eigenvectors are

Since the eigenvalues *λ*^{(1)} and *λ*^{(2)} are real and the corresponding left eigenvectors are linearly independent, the system (8) and (9) is hyperbolic and can be rewritten in characteristic form. To obtain this characteristic form, we next take the sum of (8) and (9) to obtain

and

for *n* = 1 and *n* = 2, respectively. Since (17) is identical to (9), we conclude that (9) is already in characteristic form. We now write (16) and (17) in the form

Equations (18) and (19) constitute the characteristic form of the original system (8) and (9). An advantage of (18) and (19) is that, along each family of characteristic curves, the partial differential equations have been reduced to ordinary differential equations. In regions of subsidence (i.e., where *α* = 0), information on *υ* is carried along characteristics given by (*dy*/*dt*) = 2*υ*, while information on a combination of *u* and *υ* is carried along characteristics given by (*dy*/*dt*) = *υ*. Thus, in regions of subsidence there are two distinct families of characteristics. In contrast, for regions of boundary layer pumping (i.e., where *α* = 1), the two families of characteristics become identical.

While the forcing terms *F*_{1} and *F*_{2} are too complicated to allow analytical solution of (18) and (19), the numerical solution of these ordinary differential equations can serve as the basis of the shock-capturing methods described by LeVeque (2002). In sections 5 and 6 we adopt the simpler approach of solving (1)–(7) using standard finite differences with the inclusion of horizontal diffusion to control the solution near shocks. Although this approach has some disadvantages (e.g., possible unphysical oscillations near a shock), it provides a useful guide to the expected results when global numerical weather prediction (NWP) and climate models can be run at the 100-m horizontal resolution used here.

In passing, we note that there is a less formal, more intuitive route from (8) and (9) to the characteristic forms (18) and (19). This intuitive route results from simply noting that (9) is already in characteristic form and can be directly written as (19), while the characteristic form (18) can be simply obtained by combining (8) and (9) in such a way as to eliminate terms containing the factor (1 − *α*)(∂*υ*/∂*y*).

The derivation of the characteristic form in (18) and (19) is somewhat complicated by the presence of the *w* terms in (1) and (2). Many slab boundary layer models (e.g., Shapiro 1983) neglect these *w* terms, in which case the hyperbolic nature of the problem is immediately obvious since (1) and (2) can then be written (neglecting the *K* terms) as

Even with this simplification, the three coupled, nonlinear, ordinary differential equations in (20) are too difficult to solve analytically.

## 4. Classical Ekman theory

The slab version of classical Ekman theory (Ekman 1905) is a simplification of (1)–(5). This version is obtained by neglecting the local time derivative terms, the meridional advection terms, the vertical velocity terms, and the horizontal diffusion terms. We find that the vertical velocity terms and the horizontal diffusion terms tend to be relatively small, except in regions of sharp gradients. Therefore, when we compare our solutions using the full, nonlinear slab boundary layer equations [(1)–(7)] with those produced from the classical Ekman theory equations, we can gain insight into the role of horizontal advection. With these assumptions, (1) and (2) reduce to

where the zonal geostrophic velocity *u*_{g}(*y*) is defined in terms of the specified meridional pressure gradient by

Solving the algebraic equation (21) for *u* and *υ*, we obtain

Note that (23) and (24) are actually transcendental relations for *u* and *υ* because of the dependence of *c*_{D}*U* on *u* and *υ* through (4) and (5). These transcendental equations can be solved using a variety of iterative methods.

The classical slab Ekman layer solutions [(23) and (24)] should be considered “local” in the sense that *u*(*y*) and *υ*(*y*), at a particular point *y*, depend only on the meridional pressure gradient at that point. For our present discussion, these local solutions should also be regarded as deficient because the *υ*(∂*u*/∂*y*) and *υ*(∂*υ*/∂*y*) terms have been neglected, thereby yielding solutions that are too smooth in regions of boundary layer pumping and too narrow in regions of boundary layer suction, as will be seen in the next section. While deficient, these classical Ekman solutions serve as a useful basis for comparison with the more general, nonlocal solutions presented in the next section.

## 5. Experimental results

We now specify meridional profiles of *p* and *u*_{g} for three numerical experiments. The first specified pressure field is designed to illustrate how the low-latitude, boundary layer flow responds differently to prescribed easterly versus westerly geostrophic flow along the equator. The specified pressure field is

where *b*, , and *p*_{∞} are specified constants. Substituting (25) into (22), the zonal geostrophic flow *u*_{g}(*y*) is found to be

Next, we present a second specified pressure field, which is designed to illustrate how the low-latitude, boundary layer flow responds to a more complicated prescribed geostrophic wind field. We refer to this case as the Rossby gyre case because there are two regions of easterly geostrophic flow surrounding a region of equatorial westerly geostrophic flow or westerly wind burst. The specified pressure field for the Rossby gyre case is

where *b*, , and *p*_{∞} are the same specified constants as those in (25) and (26). Substituting (27) into (22), the zonal geostrophic flow *u*_{g}(*y*) is found to be

We illustrate (25)–(28) for three experiments in Fig. 2 for the choices *p*_{∞} = 1010 hPa, *ρ* = 1.22 kg m^{−3}, *β* = 2.289 × 10^{−11} m^{−1} s^{−1}, and *b* = 1000 km. The experiments are as follows: the case of easterly geostrophic flow uses = −10 m s^{−1}, the case of westerly geostrophic flow uses = 10 m s^{−1}, and the Rossby gyre case uses = 10 m s^{−1}, resulting in the pressure perturbation hPa. Notice how the case of easterly geostrophic flow corresponds to low pressure along the equator, westerly geostrophic flow to high pressure along the equator, and the Rossby gyre case to low pressure at *y* ≈ ±707 km.

### a. Easterly and westerly geostrophic flow cases

In this subsection, we present numerical solutions of the problem (1)–(7) with the specified pressure field (25), or equivalently, the specified zonal geostrophic flow field (26). This forcing has been designed to illustrate how the nonlinear, low-latitude, boundary layer flow responds differently to prescribed easterly and westerly geostrophic flows along the equator. The numerical model uses centered, second-order spatial finite difference methods on the domain −5000 ≤ *y* ≤ 5000 km with a uniform grid spacing of 100 m and a fourth-order Runge–Kutta time differencing scheme with a time step of 5 s. The constants have been chosen as *h* = 500 m and *K* = 500 m^{2} s^{−1} for horizontal diffusivity.

In making this choice for *h*, we have taken the average value of the damping time scale observed by Deser (1993) for the boundary layer zonal flow (~7 h) and meridional flow (~17 h). For *h* = 500 m, this average damping time scale (~12 h) corresponds through (5) to the surface wind speed *U* ≈ 10 m s^{−1}. The value *h* = 500 m is also consistent with observations of the well-mixed subcloud boundary layer over the western Pacific warm pool (Johnson et al. 2001) and the eastern Pacific (Bond 1992; Yin and Albrecht 2000; McGauley et al. 2004). As these aforementioned studies discuss, the boundary layer depth varies as a function of a number of factors, especially convective activity. We have used both smaller and larger values of *h* in our model simulations to represent active and inactive convective periods in the range of 250 ≤ *h* ≤ 2000 m, and we find the results are in qualitative agreement with those using *h* = 500 m.

For the first experiment, the specified geostrophic zonal flow *u*_{g}(*y*) is centered on the equator and has an easterly maximum of 10 m s^{−1}, as shown by the dashed black curve in the left panel of Fig. 3. For the second experiment, the forcing is identical except the sign of the specified *u*_{g}(*y*) is reversed; that is, *u*_{g}(*y*) has a 10 m s^{−1} westerly maximum along the equator, as shown by the dashed black curve in the left panel in Fig. 4.

We set *u*(*y*, 0) = *u*_{g}(*y*) and *υ*(*y*, 0) = 0 as the initial conditions so, as the zonal flow slows down as a result of surface drag, the *υ*(*y*, *t*) field develops through the −*βy*(*u* − *u*_{g}) term in (2). For both of these specified *u*_{g}(*y*) forcing functions, the numerical model was integrated until a steady state in the *u* and *υ* fields was obtained, generally within 5 days. The main reason why *u* and *υ* take this long to reach a steady state is that it takes a couple of days for *υ* to strengthen since it is initially zero. The blue curves in Figs. 3 and 4 show the steady-state solutions (*t* = 120 h) for *u*(*y*, *t*), *υ*(*y*, *t*), *η*(*y*, *t*), and *w*(*y*, *t*), where the absolute vorticity *η* is related to the relative vorticity *ζ* by *η*(*y*, *t*) = *β **y* + *ζ*(*y*, *t*). We also display the solutions of the classical Ekman theory equations using (23), (24), and (26) in the red curves of Figs. 3 and 4.

For both cases, meridional flows of 2–3 m s^{−1} develop, along with large boundary layer pumping in the numerical solutions for westerly geostrophic flow (*w*_{max} ≈ 7.3 mm s^{−1} in Fig. 4) and extremely large boundary layer pumping in the numerical solutions for easterly geostrophic flow (*w*_{max} ≈ 3.2 m s^{−1} in Fig. 3). The solutions computed from classical Ekman theory produce slightly larger meridional flows with much weaker and broader boundary layer pumping (*w*_{max} ≈ 21.2 mm s^{−1} in Fig. 3 and *w*_{max} ≈ 3.2 mm s^{−1} in Fig. 4) than their numerical counterparts. Also, the classical Ekman theory solutions have larger meridional gradients of the zonal flow near the equator even though the zonal flows are of similar magnitude for both easterly and westerly cases as a result of the *β*^{2 }*y*^{2} term in the numerator of (23) requiring *u* = 0 at *y* = 0. For westerly geostrophic flow, there are regions just off of the equator that are inertially unstable (*βyη* < 0), but the numerical solutions never become inertially unstable, as shown in Fig. 4. The boundary layer suction region in the classical Ekman theory solutions is too narrow and is too intense (*w*_{min} ≈ −21.2 mm s^{−1} in Fig. 4) along the equator in the westerly case. In fact, the classical Ekman solutions of *u*, *υ*, *w* for the easterly and westerly geostrophic flow cases are exactly the same except they are of opposite signs. These discrepancies between the numerical model and classical Ekman theory solutions are mainly attributed to the *υ*(∂*υ*/∂*y*) term in the meridional equation of motion. This term acts to sharpen and strengthen regions of boundary layer pumping and broaden and weaken regions of boundary layer suction.

We now interpret the easterly and westerly numerical model experiments in terms of the simplified dynamics [(20)], written in the form

where

is the absolute angular momentum and (*d*/*dt*) = (∂/∂*t*) + *υ*(∂/∂*y*). Consider first the numerical solutions for the case of easterly geostrophic flow, which is somewhat more straightforward to interpret dynamically and is shown in the blue curves of the four panels of Fig. 3. In this case, the −(*c*_{D}*Uu*/*h*) term in the zonal equation of motion immediately begins to slow down the boundary layer easterly flow at a rate of approximately 14.6 m s^{−1} day^{−1}. Because of this slowing of the zonal flow, the −*βy*(*u* − *u*_{g}) term in (30) becomes negative for *y* > 0 and positive for *y* < 0. This subsequently produces a 3.1 m s^{−1} northerly flow just north of the equator and a 3.1 m s^{−1} southerly flow just south of the equator. These opposing meridional flows tighten up until a near discontinuity, or shocklike structure, in *υ* is produced at the equator. Associated with this shocklike structure in *υ* is a near singularity in the boundary layer pumping *w*. There is no singularity in the absolute vorticity *η*, although there is a near discontinuity in *η* at the equator associated with the kink in the *u* field. For this case of easterly geostrophic flow along the equator, horizontal diffusion is the main physical process balancing the *υ*(∂*υ*/∂*y*) term using our current numerical methods (see discussion of Fig. 7); therefore, we refer to this narrow boundary layer pumping region as shocklike.

The numerical solutions for the case of westerly geostrophic flow are shown by the blue curves in the four panels of Fig. 4. In this case, there is high pressure along the equator, so that, as zonal surface friction begins to slow the boundary layer westerly flow, the −*βy*(*u* − *u*_{g}) term in (30) produces a 2.9 m s^{−1} meridional flow that is divergent near the equator. There are two narrow regions of enhanced boundary layer pumping at *y* ≈ ±950 km, a double ITCZ (Zhang 2001). Two small regions of slightly supergeostrophic boundary layer zonal flow are produced in the regions 875 < |*y*| < 975 km. To better understand this feature, we analyze the effect of the meridional flow on *u* through the relation

which is obtained from the definition of *m*. The *βyυ* term in (32) is positive everywhere since *υ* > 0 where *y* > 0 and *υ* < 0 where *y* < 0. Therefore, for large enough *βyυ*, we obtain from (32) the result that (*du*/*dt*) > 0 even though the surface friction term in (29) requires (*dm*/*dt*) < 0. The solution for *υ* resembles the N-wave solution found for Burgers’s equation (Whitham 1974), but the magnitude of the maximum *w* (at *y* ≈ ±950 km) is much smaller than the maximum of *w* in the easterly wind case. A true Burgers’s equation N-wave structure has been arrested by the combination of the −(*c*_{D}*Uυ*/*h*) and −*βy*(*u* − *u*_{g}) terms in (30). Unlike the case of easterly geostrophic flow there are two narrow regions of anomalous absolute vorticity *η* located at nearly the same locations as the two peaks in enhanced pumping. The maxima in *η* are displaced slightly poleward of the peak in *w*, near *y* ≈ ±1000 km. This displacement between *ζ* and *w* is in general agreement with solutions for the tropical cyclone boundary layer (Williams et al. 2013).

As discussed in section 3, characteristics are very helpful in understanding the development of shocks. Since the slab boundary layer model used here contains horizontal diffusion, we do not strictly have a hyperbolic system. However, for the cases discussed here, we can consider the system as effectively hyperbolic since the diffusion terms are either negligible everywhere or nonnegligible only in a small region (see Figs. 7–9, 12, and 13). As can be seen in (18) and (19), the two families of characteristics are given by (*dy*/*dt*) = *υ* and (*dy*/*dt*) = (2 − *α*)*υ*, so that in regions where *w* > 0 (*α* = 1) the two families are identical and are the same as particle trajectories. In regions where *w* < 0 (*α* = 0), the two families of characteristics are distinct, and one family is the same as particle trajectories while the other family moves twice as fast as particle trajectories in the meridional direction. In the following figures, we present particle trajectories, which can be interpreted as both families of characteristics in regions where *w* > 0 and as one of the two families in regions where *w* < 0.

Figures 5 and 6 illustrate trajectory lines along with filled contours of *u*(*y*, *t*), *υ*(*y*, *t*), and *w*(*y*, *t*) for the cases of easterly and westerly geostrophic flow, respectively. The trajectories were computed by numerically integrating (*dy*/*dt*) = *υ* using the same 5-s time step used for the numerical solutions of (1) and (2). In the right panel of Figs. 5 and 6, we zoom in on the region of sharp gradients (−500 ≤ *y* ≤ 500 km) shown in the left two panels of Fig. 5 and (500 ≤ *y* ≤ 1500 km) in the left two panels of Fig. 6. Also, recall from (18) and (19) that regions of boundary layer pumping contain only one family of characteristics while regions of boundary layer suction contain two families of characteristics. Therefore, the trajectories we show are essentially the same as characteristics in regions where *w* > 0.

Initially, the particles are stationary, owing to the initial condition *υ* = 0. As time progresses, the trajectories begin to turn toward the low pressure region along the equator in Fig. 5 and away from the high pressure region along the equator in Fig. 6. In Fig. 5, the trajectories intersect at the equator, producing a shocklike structure at *t* ≈ 55 h. In contrast, the trajectories in Fig. 6 converge, but they never intersect, suggesting the absence of any shocklike structures. Instead, two narrow regions of enhanced boundary layer pumping form at *y* ≈ ±950 km and *t* ≈ 55 h. Both cases take approximately 2 days to form narrow regions of boundary layer pumping because it takes some time for the meridional winds *υ* to build up in magnitude. The regions of largest boundary layer subsidence, which are away from the equator in the easterly flow case and along the equator in the westerly flow case, do not take as long as the regions of boundary layer pumping to become established.

During the first 40 h of the easterly flow simulation, the left panel of Fig. 5 illustrates that (*du*/*dt*) > 0 (slowing down of the easterlies) because of the zonal surface drag, but after that time (*du*/*dt*) < 0 (easterlies speed up) because of the increase in the magnitude of *βyυ*. Despite this, supergeostrophic flow is never produced in the easterly flow case, unlike the westerly flow case. The left panel of Fig. 6 illustrates that *u*(*y*, *t*) decreases because of the presence of zonal surface drag, until about 50 h. At this time, *u*(*y*, *t*) begins to increase along trajectories despite the continual decrease in *m*(*y*, *t*). This increase in *u*(*y*, *t*) to supergeostrophic speeds (*u* > *u*_{g}) is associated with the decrease in *υ*(*y*, *t*) and shocklike structures (see Fig. 9).

To further assist in our comparison between the cases of easterly and westerly geostrophic zonal flow, we compute the steady-state values of individual terms in (1) and (2). First, we rewrite (1) and (2) in the form

where we have used *η* = *βy* − (∂*u*/∂*y*) and (22). We illustrate the individual terms for the easterly and westerly geostrophic flow cases in Figs. 7–9. In the left panels of Figs. 7 and 9, the terms on the right-hand side of (33) are shown, and the right panels of Figs. 7 and 9 show the terms on the right-hand side of (34). Figure 8 zooms in on Fig. 7 in the region of sharp gradients, −5 ≤ *y* ≤ 5 km. Note that the blue and black curves in the right panels of Figs. 7 and 9 are the same in regions of subsidence since *w*(1 − *α*)*υ*/*h* = −*υ*(∂*υ*/∂*y*).

Figure 7 shows that both horizontal diffusion terms are negligible away from the equator but are very large near the equator, especially in the meridional momentum equation, where *K*(∂^{2}*υ*/∂*y*^{2}) ≈ ±500 m s^{−1} day^{−1}, as shown in Fig. 8. In this near-equatorial region, *K*(∂^{2}*υ*/∂*y*^{2}) balances −*υ*(∂*υ*/∂*y*), signifying the presence of a shocklike structure. In the left panel of Fig. 7, *K*(∂^{2}*u*/∂*y*^{2}) balances −(*c*_{D}*Uu*/*h*) near the equator since *υη* → 0.

The left and right panels of Fig. 9 illustrate the lack of horizontal diffusion, suggesting the absence of shocklike structures. More specifically, the right panel of Fig. 9 illustrates that where −*υ*(∂*υ*/∂*y*) is largest (*y* ≈ ±900 km), it is balanced by the −(*c*_{D}*Uυ*/*h*) and −*βy*(*u* − *u*_{g}) terms, as suggested previously. In this way, a shock or shocklike structure is avoided. Even without shocklike structures, the classical Ekman theory solutions shown in Fig. 4 cannot accurately reproduce the numerical solutions. This is mainly due to the lack of the *υ*(∂*υ*/∂*y*) term in classical Ekman theory that plays an important role in sharpening and strengthening boundary layer pumping regions and broadening and weakening boundary layer suction regions.

### b. Rossby gyre case

In the third experiment, we compute numerical solutions of the problem (1)–(7) with the specified pressure field from (27), or equivalently, a geostrophic zonal flow field with easterly winds poleward of a westerly wind burst along the equator, as given in (28). We refer to this as the Rossby gyre case because it resembles the meridional structure of both theoretical and observed Rossby waves (Matsuno 1966; Kiladis et al. 1994; Kiladis and Wheeler 1995). The geostrophic zonal winds *u*_{g}(*y*) are centered on the equator with a westerly maximum of 10 m s^{−1}, as shown in the dashed black curve in the left panel of Fig. 10. The blue curves in Fig. 10 show the steady-state solutions (*t* = 120 h) for *u*(*y*, *t*), *υ*(*y*, *t*), *η*(*y*, *t*), and *w*(*y*, *t*). We also display the solutions of the classical Ekman theory using (23), (24), and (28) in the red curves of Fig. 10.

Meridional flows of 2 m s^{−1} develop, along with large boundary layer pumping in the numerical solutions (*w*_{max} ≈ 26 mm s^{−1} in Fig. 10). The solutions computed from classical Ekman theory produce slightly larger meridional flows with much weaker and broader boundary layer pumping (*w*_{max} ≈ 5.8 mm s^{−1} in Fig. 10) than the numerical model. Also, the classical Ekman theory solutions have larger meridional gradients of the zonal flow near the equator and produce inertially unstable (*βyη* < 0) regions near the equator. Once again, the differences between the numerical model and classical Ekman theory solutions are mainly attributed to the *υ*(∂*υ*/∂*y*) term in the meridional equation of motion.

Similar to the case of only westerly geostrophic flow, there is low pressure away from the equator; therefore, the −*βy*(*u* − *u*_{g}) term in (30) produces a meridional flow of approximately 2 m s^{−1} that is divergent near the equator. Unlike the westerly geostrophic flow case, there are pressure minima centered at *y* ≈ ±707 km; therefore, the −*βy*(*u* − *u*_{g}) term in (30) produces equatorward flow of approximately 1.2 m s^{−1} poleward of *y* ≈ ±707 km. There are two narrow regions of boundary layer pumping in either hemisphere, a double ITCZ, similar to the westerly geostrophic flow case. The two narrow regions of boundary layer pumping are stronger and narrower in this case; peak values are over 3 times as large as the case of only westerly geostrophic flow along the equator. The main reason why the pumping is stronger and more concentrated in the Rossby gyre case is that the region of westerly geostrophic flow along the equator is narrower than in the westerly geostrophic flow case. Because of this, the regions of maximum rising motion are located closer to the equator (*y* ≈ ±620 km) compared to the westerly geostrophic flow case (*y* ≈ ±950 km). The magnitude of supergeostrophic boundary layer zonal winds is larger in the Rossby gyre case than the westerly geostrophic flow case. Also, the second derivative of the steady-state westerly flow along the equator is negative (*d*^{2}*u*/*dy*^{2} < 0) as opposed to positive (*d*^{2}*u*/*dy*^{2} > 0) in the westerly case. Along with two peaks in boundary layer pumping, there are two narrow regions of enhanced absolute vorticity *η* displaced slightly poleward of the peak in boundary layer pumping, at *y* ≈ ±625 km.

We illustrate trajectory lines along with filled contours of *u*(*y*, *t*), *υ*(*y*, *t*), and *w*(*y*, *t*) in Fig. 11 for the Rossby gyre case. In the right panel in Fig. 11, we zoom in on the region of sharp gradients (0 ≤ *y* ≤ 1000 km) shown in the left panels of Fig. 11. The behavior of this case is very similar to the case of only westerly geostrophic flow with two regions of enhanced pumping forming in both hemispheres but around *y* ≈ ±620 km instead of *y* ≈ ±950 km. The trajectories associated with the westerly geostrophic flow do converge closer in the Rossby gyre than the case of westerly geostrophic flow, but once again they never intersect. The two narrow regions of enhanced boundary layer pumping form after the sufficient increase in the meridional winds at *t* ≈ 55 h. The trajectories poleward of the regions of boundary layer pumping, where the zonal flow is easterly, turn equatorward and converge as well, but they do not intersect. The *w*(*y*, *t*) field has finer-scale structures poleward of the maxima in boundary layer pumping in the Rossby gyre case than the other two cases.

We illustrate the steady-state values of individual terms in (33) and (34) for the Rossby gyre case in Fig. 12 and Fig. 13. The left panel of Fig. 12 shows the terms on the right-hand side of (33); the right panel of Fig. 12 shows the terms on the right-hand side of (34). Figure 13 zooms in on the Northern Hemisphere region of sharp gradients, 612 ≤ *y* ≤ 622 km, in the right panel of Fig. 12. Recall that the blue and black curves in the right panel of Fig. 12 are the same in regions of subsidence since *w*(1 − *α*)*υ*/*h* = −*υ*(∂*υ*/∂*y*).

Both panels of Fig. 12 illustrate that horizontal diffusion is negligible everywhere except close to the two narrow peaks in boundary layer pumping (*y* ≈ ±620 km). However, Fig. 13 illustrates that where −*υ*(∂*υ*/∂*y*) is largest, it is balanced mainly by the −(*c*_{D}*Uυ*/*h*) and −*βy*(*u* − *u*_{g}) terms and not horizontal diffusion. We do not define this feature as a shocklike structure because even though horizontal diffusion is playing an important role near the sharp gradients, it is not the physical process that prevents a shock. Although we believe that horizontal diffusion is not vital to preventing a shock, spurious oscillations form in the time evolution of the wind fields if horizontal diffusion is discarded. This examination of the individual terms provides a better understanding of the roles that horizontal diffusion, surface friction, and ageostrophic zonal flow play in forming narrow ITCZs.

## 6. Concluding remarks

The boundary layer wind field near the ITCZ has been interpreted in terms of a zonally symmetric slab boundary layer model. The narrowness of the ITCZ has been explained by the formation of narrow regions of enhanced boundary layer pumping, associated with the dynamical role of the *υ*(∂*υ*/∂*y*) term in the meridional momentum equation. In the case of easterly geostrophic flow (low pressure along the equator), a shocklike structure is produced at the equator; that is, the meridional flow *υ* would become discontinuous and the boundary layer pumping would become singular (*w* → ∞ at *y* = 0) if *K* = 0. In the two cases of westerly geostrophic flow (high pressure along the equator), true shocks are not produced, but narrow regions of enhanced boundary layer pumping are produced on both sides of the equator.

When the numerical solutions are compared to those of classical Ekman theory, the important role of the horizontal advective terms *υ*(∂*u*/∂*y*) and *υ*(∂*υ*/∂*y*) becomes evident. The solutions of the classical Ekman theory tend to be too smooth in regions of boundary layer pumping, especially near the equator. They also have deficiencies in producing regions of boundary layer suction that are both too strong and concentrated when compared to the numerical solutions. Although the numerical solutions from the slab boundary model are based on the assumption of zonal symmetry, they may have relevance to the boundary layer of the Madden–Julian oscillation (MJO). This is because the horizontal winds surrounding the MJO convective envelope are westerly to the west and easterly to the east (Schubert and Masarik 2006). We speculate that there may be two regions of enhanced boundary layer pumping almost symmetric about the equator on the west side and one on-equatorial convective region on the east side of the MJO convective envelope. We anticipate that the two regions of enhanced boundary layer pumping on the west side are not as narrow as the convective region on the east side. However, the boundary layer and free-tropospheric thermodynamics may respond to easterlies along the equator by producing convection near the equator before the convergence becomes as thin as the steady-state solutions presented in the slab model suggest, especially because a narrow ITCZ centered on the equator is rarely observed in nature. It is also possible that the assumption of zonal symmetry breaks down for a phenomenon such as the MJO. However, narrow or shocklike structures associated with zonal and meridional advection may both become important in this scenario.

The interactions between boundary layer convergence and free-tropospheric convection in and near the ITCZ involve complexities that our slab model cannot address. However, there are numerous studies that analyze the interactions between low-level dynamics and thermodynamics in and near the ITCZ, such as Lindzen and Nigam (1987), Waliser and Somerville (1994), Tomas and Webster (1997), Liu and Moncrieff (2004), Gu et al. (2005), Raymond et al. (2006), and Sobel and Neelin (2006). Waliser and Somerville (1994), Liu and Moncrieff (2004), and Sobel and Neelin (2006) place emphasis on processes involving nonlinear dynamics and thermodynamics, which may be most relevant near the equator where classical Ekman theory breaks down. On-equatorial atmospheric boundary layer convergence may not typically couple to convection as a result of the upwelling of cold ocean water (Charney 1969, 1971; Pike 1971, 1972; Mitchell and Wallace 1992; Liu and Xie 2002), particularly under easterly winds. Also, regions of upwelling are typically associated with large boundary layer static stability (Bond 1992; Yin and Albrecht 2000) that can limit the efficiency of boundary layer pumping (Gonzalez and Mora Rojas 2014). Nonetheless, the high-resolution slab model simulations of single and double ITCZs are intriguing, especially as global models continue to resolve finer-scale features.

For the boundary layer structures simulated here, we choose to use the term “boundary layer shock” and reserve the term “front” for structures that arise not from *υ*(∂*υ*/∂*y*) but rather from the combination of *υ*(∂*u*/∂*y*), *w*(∂*u*/∂*z*), *υ*(∂*θ*/∂*y*), and *w*(∂*θ*/∂*z*), with the rotational flow *u* and the potential temperature *θ* being related by thermal wind balance. This terminology helps distinguish features that can be accurately modeled using the geostrophic balance assumption (i.e., fronts) from features that cannot be modeled using geostrophic balance (i.e., boundary layer shocks).

The phenomena of narrow and shocklike ITCZs put demanding horizontal resolution requirements on global NWP and climate models that are as strict as those for accurate simulation of moist convection. In view of the importance of boundary layer pumping in determining the location of diabatic heating, accurate ITCZ simulations probably require accurate simulations of such finescale aspects of the boundary layer. In closing, we reiterate that the ITCZ boundary layer is an environment conducive to the formation of zonally elongated regions of enhanced boundary layer pumping. These regions are an important factor in determining the organization of deep convection. However, enhanced boundary layer pumping due to surface friction does not completely dictate the location of deep convection since other factors, such as a dry midtropospheric environment, can lead to convective downdrafts that produce expanding cold pools and hence also help determine the organization of convection (Khairoutdinov et al. 2010). Considerable work remains to understand the coupling of these complex processes.

## Acknowledgments

We thank Paul Ciesielski, Thomas Birner, Eric Maloney, and Donald Estep for their insightful discussions. We would like to acknowledge the reviewers, Adam Sobel and David Raymond, for their constructive comments on the manuscript. This research has been supported by the National Science Foundation under Grant AGS-1250966 and under the Science and Technology Center for Multi-Scale Modeling of Atmospheric Processes, managed by Colorado State University through Cooperative Agreement ATM-0425247.

### APPENDIX

#### Derivation of the Boundary Layer Equations

The starting point in the derivation of (1) is the conservation relation for the absolute angular momentum. On the equatorial *β* plane, the absolute angular momentum is *m* = *a*[*u* + (1/2)*β*(*a*^{2} − *y*^{2})], which is an approximation of the spherical version *a*[*u* cos*ϕ* + (1/2)*β*(*a*^{2} − *a*^{2} sin^{2}*ϕ*)]. This conservation relation can be written in the flux form

where *m*_{g} = *a*[*u*_{g} + (1/2)*β*(*a*^{2} − *y*^{2})] is the absolute angular momentum of the geostrophic flow above the boundary layer. According to (A1), there are five processes that can cause changes in the boundary layer absolute angular momentum: (i) meridional divergence of the meridional advective flux, (ii) upward flux of *m* when *w* > 0, (iii) downward flux of *m*_{g} when *w* < 0, (iv) loss of *m* through surface drag, and (v) meridional divergence of the meridional diffusive flux. To convert (A1) into a more convenient form we differentiate the second term as a product and then make use of the continuity equation [(3)] to obtain (1).

The derivation of the meridional momentum equation [(2)] proceeds in a similar fashion. The flux form is

The pressure gradient force in the boundary layer is now assumed to be equal to the pressure gradient force in the region above the boundary layer, where geostrophic balance exists. To convert (A2) into a more convenient form we differentiate the second term as a product and then make use of the continuity equation [(3)] to obtain (2).

## REFERENCES

*Proc. WMO/IUGG Symp. on Numerical Weather Prediction*, Tokyo, Japan, Japanese Meteorological Agency,

*Mathematical Problems in the Geophysical Sciences*, W. H. Reid, Ed., Lectures in Applied Mathematics, Vol. 13, American Mathematical Society, 355–368.

*J. Adv. Model. Earth Syst.*,

**1**, 15, doi:.

*Finite Volume Methods for Hyperbolic Problems.*Cambridge University Press, 558 pp.

*Linear and Nonlinear Waves.*John Wiley and Sons, 363 pp.

*J. Meteor. Soc. Japan*,

**49A**, 691–698.

## Footnotes

Current affiliation: Joint Institute for Regional Earth System Science and Engineering, University of California, Los Angeles, Los Angeles, and Jet Propulsion Laboratory, California Institute of Technology, Pasadena, California.