The Dynamical Role of the Chesapeake Bay on the Local Ozone Pollution Using Mesoscale Modeling—A Case Study

: This study investigated the dynamic inﬂuence of the Chesapeake Bay (CB) on local ozone (O 3 ) concentration and distribution using a weather forecasting model. The Weather Research and Forecasting model coupled with Chemistry (WRF–Chem) was employed to simulate O 3 production and transportation near the CB. Baseline (water) as well as sensitivity (nowater) model experiments of bay circulation were carried out with and without bay water by changing the water surface from water to land (loam). First, the model performance simulating O 3 was evaluated using the baseline experiment results and AirNow surface wind and O 3 observations. The results showed that the model overestimates surface O 3 by up to 20–30%. Further, the comparisons of the baseline and sensitivity experiments revealed higher O 3 mixing ratios, primarily due to the resulting bay breeze circulation. These increases, after considering model overestimation, represent a mean bay dynamics circulation-induced contribution of up to 10% at night and 5% during the day. Furthermore, the boundary layer over northern CB, where it is at its narrowest width, was higher (by 1.2 km on average) during daytime due to higher surface temperatures observed. The boundary layer depth difference between the northern, central, and southern regions of the bay leads to a differential in the role of bay circulation dynamics in the observed O 3 increase. The relatively wider swath of water surface over southern CB resulted in a lower boundary layer depth and stronger breeze circulation and this circulation contributed to O 3 concentrations. Moreover, since the case selected has a minimal bay breeze circulation, the associated surface ozone enhancements represent what is expected at least at a minimum.


Introduction
Tropospheric (ground-level or surface) ozone (O 3 ) has significant environmental and human health impacts. Long-term exposure to high O 3 concentration air leads to serious health issues, such as irritated lungs, aggravated bronchitis, emphysema, and asthma [1][2][3]. O 3 also has significant impacts on crop yields, forest growth, and species composition [4]. Recent studies have focused on quantifying the impact of current O 3 concentration on net ecosystem production as compared to the conditions of the pre-industrial era. [5,6]. To mitigate surface pollution in the U.S., the Clean Air Act requires the Environmental Protection Agency (EPA) to set National Ambient Air Quality Standards (NAAQS) [7]. The current NAAQS, in effect since 2015, sets an ambient 8 h O 3 mixing ratio criterion of 70 parts per billion by volume (ppbv) [8,9]. In particular, the O 3 mixing ratio over Atmosphere 2022, 13, 641 2 of 20 some coastal regions exceeds the NAAQS more frequently compared to other inland regions. Examples of such regions include the Chesapeake Bay (hereafter referred to as CB) [10,11], Galveston Bay [12], the Great Lakes [13][14][15][16], and the Great Salt Lake [17]. All these locations are in close proximity to water bodies surrounded by land, and related research is ongoing to investigate the causes of the exceedance of O 3 concentration [11,18,19]. As a secondary pollutant, ground-level O 3 depends on its precursors, chemical reactions, and meteorological conditions [20][21][22]. Understanding them helps to track O 3 pollution sources and sinks and is critical for improving the accuracy of O 3 pollution forecasting.
Field campaigns provide insights into physical and chemical processes in the atmosphere under a diverse set of meteorological conditions over the aforementioned locations. A number of campaigns have been conducted to investigate bay/lake effects on local meteorology and air quality, e.g., Deriving Information on Surface Conditions from Column and Vertically Resolved Observations Relevant to Air Quality from 27 June to 31 July 2011 (DISCOVER-AQ) [23,24], the Border Air Quality and Meteorology Study from 20 June to 10 July 2007 (BAQS-Met) [16], the Ozone Water-Land Environmental Transition Study from 5 July to 3 August 2017 (OWLETS) and from 6 June to 6 July 2018 (OWLETS2) [11,18], and the Lake Michigan Ozone Study (LMOS) [13,14,19,25]. Findings on the impacts of CB on O 3 from both dynamical and chemistry perspectives have been generated and a link to the bay/lake water body has also been suggested by these field campaign studies [10,13,14,17,19,[23][24][25].
One of the most important contributions made by the previous studies concerns the role of dynamics originating from the local bay/lake and coastal circulations. The bay/lake dynamical effects contributing to the local meteorology and air quality are: (1) the developing land/bay breeze and resulting diurnal dynamics that alter the near-shore meteorology through small-scale circulations, e.g., the bay breeze yields high pressure over CB and low pressure inland near CB [26]; (2) associated changes in cloud formation, which, in turn, modify the radiation environment, thereby causing shielding for O 3 photolysis [10], e.g., the presence of clear skies over CB and later development of cumulus clouds due to higher surface temperatures and elevated boundary layers [26]; (3) changes in the vertical mixing of gases, associated with changes in plume updraft velocity over land compared to water, which impact O 3 concentration, e.g., the development of a bay breeze during the late morning and early afternoon along the western shore of CB was found to lead to a convergence of pollutants [10]; and (4) modification of the inversion layer (strong inversion over land at night as opposed to over the bay/lake) and its impact on the dilution and venting of gases.
Most of these previous studies utilized observations from field campaigns to analyze the effects of the bay/lake and evaluate model performance by comparing model results with observations [10,25,26]. To our knowledge, not many studies have taken advantage of recent modeling improvements, especially the online coupling of meteorology and chemistry, to specifically address the dynamical role of CB on O 3 concentrations.
Even though some previous reports have discussed the numerous possible dynamic and chemical processes that alter near-shore pollution, they have focused on model resolution sensitivity and model performance evaluation [10,11,18,[25][26][27][28][29]. Chief among the dynamical processes are the development of breeze (bay and land) and the differential sensible heating resulting in the land-bay horizontal temperature gradients. The modelling of sea breezes depends on the land surface sensible heat flux, ambient geostrophic wind, atmospheric stability and moisture, water body dimensions, terrain height and slope, Coriolis parameter, surface roughness length, and shoreline curvature [30]. To capture these features in a model simulation, spatial resolution is a key factor to determine model performance when simulating different scales. Studies have demonstrated that finer horizontal spatial resolution (0.5 km [10]; 1.33 km [27]) model runs agreed better with pollutant observations at the top of the boundary layer and resolved the structure of the bay breeze better as compared to coarser model resolutions (13.5 km, [10]; 12 km, [27]). Similar results were also reported by Jiménez et al. (2006), who concluded that O 3 simulation during sea breezes is better with a fine resolution model than with a coarse resolution model [31].
This study further investigates the influence of CB on the O 3 mixing ratio distribution by employing a mesoscale model, the Weather Research and Forecasting model coupled with Chemistry (WRF-Chem) [32]. It aims to quantify the relative influence of the bay-land dynamics on the O 3 concentration near bay areas. One baseline experiment ("water") and one sensitivity experiment ("nowater") were performed with the normal configuration and by altering the model configuration of the underlying water over CB to land (loam), respectively. First, the baseline model simulation performance was evaluated by the surface O 3 mixing ratio observations to make sure that model simulations were reliable. Then, both horizontal and vertical O 3 distribution differences between "water" and "nowater" over CB were analyzed to study the resulting dynamics and their influence on the temporal spatial O 3 concentration distribution.
This paper is organized as follows: descriptions of the model configuration, study domain, background, case study, and datasets are given in Section 2, followed by a model performance evaluation in Section 3. The results and analysis are presented in Section 4, and concluding remarks are made in Section 5.

Model Configuration
WRF-Chem (version 3.7) integrates the air quality component consistent with the meteorological components [32]. The gas-phase chemistry and aerosol module were based on the Carbon Bond Mechanism Version Z (CBM-Z) [33] and the Model for Simulating Aerosol Interactions and Chemistry (MOSAIC) [34], respectively. The Yonsei University (YSU) planetary boundary layer scheme [35,36] was used here based on the model performance evaluation and recommendation for such an area [37]. Radiation treatment utilized the Rapid Radiative Transfer Model for General Circulation Model (RRTMG) short-wave and long-wave radiation schemes [38]. The radiation scheme was in line with what was recommended for the model simulations over the continental U.S. by the WRF developer team [39]. Figure 1 shows the model domain, with the outer study area set to include the eastern U.S. The outer domain was further resolved by one two-way nested grid region that progressively focused and centered on the study region, CB. The outer domain (d01) was set at 9 km × 9 km horizontal spatial resolution, while the nested inner domain (d02) had a higher horizontal spatial resolution (3 km × 3 km) because it exclusively focused on the detailed investigation of O 3 concentration and its evolution. Thirty-five vertical levels were used with about 15 levels set below 800 hPa to resolve the PBL, and the model top was at 100 hPa.

Topography of the CB Surrounding Region
CB is located on the western coast of the Atlantic Ocean and surrounded by the states of Maryland, Virginia, and Delaware. It is approximately 320 km long from north to south, 4.5 km wide at its narrowest point from west to east, and 48 km at its widest point. The average depth is 6.4 m, with a maximum of 53 m [41]. The terrain towards the east of CB is flat, and its terrain height is less than 50 m. Mountains are located 100 km west of CB, with the terrain height more than 1000 m ( Figure 1). CB alters moisture, wind speed and direction, and atmospheric stability. Due to the evaporation of the water body over CB, the air moisture level near it is higher than in the surrounding area. It has a relatively high temperature at night, which makes the atmosphere unstable, while, in the day, it has a relatively low temperature, which stabilizes the atmosphere.

Case Study
There are several points that need to be considered that may have an impact on the choice of the selection of 3 June 2015. First, O 3 is mostly a summertime problem, and early June is the beginning of summer and a heightened pollution period in this region. Stauffer and Thompson (2015) used O 3 climatology data for Baltimore, MD, and Hampton, VA, from May-September, 1981-2010, and distinguished three day types: no bay breeze, rainy/cloudy, and bay breeze [42]. The O 3 concentration ranges were 40-100 ppbv (median at 60 ppbv) on no bay breeze days, 20-80 ppbv (median at 50 ppbv) on rainy/cloudy days, and 40-100 ppbv (median at 70 ppbv) on bay breeze days. The O 3 diurnal variation varies from 15 ppbv at night to 70 ppbv in the day. The 8 h O 3 mean on bay breeze days is 3-5 ppbv higher than on no bay breeze days. Second, this provides for a relatively uncomplicated pollution day, a time just after an exit southward of a cold front (Figure 2a), providing a relatively less polluted day and the start of a bay-land interaction in the PBL development. Although this case will be slightly different from a muggy and stagnant case, it will serve as an important case for including all the physical processes that may be at work in the bay-land interaction. However, there was a low-pressure system near North Carolina which formed clouds and precipitation over the study region (Figures 2b and 3). The cloud fraction was 100% from the MODIS cloud product (Figure 3c,d). The 24 h precipitation map shows that the rainfall was around 0.1-0.5 in (or 0.25-1.25 cm, Figure 2b). This means that less solar radiation was expected over CB, which would reduce the photochemical reactions of O 3 and its precursors. So, the chemistry influence on O 3 generated was limited due to the absence of sunlight. Third, the prevailing winds on this day were northeasterly rather Atmosphere 2022, 13, 641 5 of 20 than the typical northwesterly. This would avoid the influence of pollutants transported from the Ohio River Valley. In addition, the synoptic scale and rainy/cloudy conditions on this day constituted minimal impacts of CB on meteorology and O 3 formation. This case will capture the physical processes involved and their impacts on controlling the mesoscale distribution of pollutants in this region and serve as a guide to future studies. An effort is under way to repeat this work for several cases and synoptic wind types in the future.

Datasets
Relevant information on the data used in the simulation and sensitivity studies as well as model evaluation is presented in Table 1. Table 1. Data used as inputs to WRF-Chem simulations.

Data Type Description
Meteorological initial and boundary conditions Northern American Regional Reanalysis (NARR) dataset, which is a high-resolution model-assimilated observation dataset from National Centers for Environmental Prediction (NCEP). The NARR covers the time period from 1979 to near present and provides 3-hourly and monthly data at a resolution of approximately 32 km with 29 pressure levels, from 1000 to 100 hPa. Both land use and water land mask datasets are used in this study. The land use (soil type) dataset from the United States Geological Survey (USGS) has 16 soil categories ( Figure 4). Since the spatial resolution of the land mask dataset is 30 arcsec, it cannot resolve the fine streams and branches feeding into CB. As a result, the shape of the water mask does not exactly follow the CB outline, which may lead to uncertainties about the size of CB. In the sensitivity experiment, the land mask was altered from water to land over CB (Figure 4a,b). The land use index was modified to loam, since the land over the eastern shore of CB is loam (Figure 4c,d). The CB surface temperature was from NARR reanalysis data. The sea surface temperature (SST) data are from large-scale forcing data, which are from the NARR dataset. The EPA AirNow program provides forecasts and near real-time observed air quality information across the U.S., Canada, and Mexico (http://www.epa.gov/AirNow, accessed on 20 September 2020). It receives air quality observations from over 1000 monitoring stations and collects forecasts for more than 300 cities. For this study, wind and O 3 mixing ratios measured near the surface were used as evaluation datasets to estimate model simulation performance. In the following comparison, we only considered the O 3 mixing ratio difference between baseline and sensitivity experiments, except for the model performance evaluation, where the O 3 mixing ratio was used.

Model Simulation Evaluation
To stabilize the meteorological and chemical fields, the baseline experiment was conducted from 1200 UTC, 1 June, to 0000 UTC, 4 June 2015. The first 36 h simulations were treated as spin-up and the remaining 24 h simulations (from 0000 UTC, 3 June, to 0000 UTC, 4 June) were selected for analysis. The sensitivity experiment was the same as the baseline experiment, except for the surface land type over CB. The surface land type was changed from water to land (loam) on 0000 UTC, 2 June 2015. The sensitivity experiment continued to run until 0000 UTC, 4 June 2015. In the following comparison between the baseline and sensitivity experiments, we selected the simulations from 0000 UTC, 3 June, to 0000 UTC, 4 June 2015.
The wind patterns influenced the instantaneous surface O 3 mixing ratio horizontal spatial distribution. Both model-simulated winds and AirNow observed winds show that the prevailing winds were northeasterly in the northeastern domain, easterly over mountains, northerly over Virginia state, and southerly over the Atlantic Ocean in the southeastern domain ( Figure 5). For example, the high O 3 mixing ratio (around 35 ppbv) over the east of the Appalachian Mountains was primarily due to the accumulated O 3 transported by the northeast winds from the polluted region nearby. When O 3 reached the eastern side of the mountains, it was blocked, accumulated, and lofted, leading to an apparent increase in O 3 concentration ( Figure 5). On the other hand, the uniform northeasterly wind (5-15 m s −1 ) over the relatively flat coastal terrain of Delaware-Maryland-Virginia favored a uniform surface O 3 mixing ratio spatial distribution. At the northern CB, the model-simulated O 3 mixing ratio agreed with the AirNow O 3 observation well. The relatively high O 3 concentration (>25 ppbv) existed over the south of CB. The high wind speed transported O 3 from CB, and local stagnant wind and weak convergence also contributed to the high O 3 concentration. Even though the influence of prevailing winds on the O 3 concentration distribution was important, the sensitivity experiment was able to capture the O 3 concentration differences due to the dynamics of CB. This is different from the study conducted by Stauffer et al. (2015), in which conclusions were based on the condition that the mesoscale or synoptic-scale wind must be absent [41]. The model-simulated surface winds and AirNow surface winds agree in both speed and directions over the study domain. In comparison to the AirNow surface O 3 concentration, the model results overestimated the surface O 3 mixing ratio over the middle and west of the study region, while they underestimated O 3 in the Delaware and New Jersey regions due to the clean air from the Atlantic Ocean. Overall, the model simulation showed O 3 mixing ratio surface spatial distribution patterns that agreed well with observations ( Figure 5). We investigated the hourly surface O 3 mixing ratio diurnal evolution for the model and observations to evaluate the model performance on capturing the O 3 diurnal variation. To compare the model simulation with the observations, the simulated O 3 mixing ratios in nine 3 km × 3 km grids surrounding the AirNow station were aggregated to one 9 km × 9 km grid region and their means and standard deviations were calculated. Figure 6 shows hourly surface O 3 mixing ratio diurnal variations from three AirNow stations (Millington, Fair Hill, and Aldino) close to the north of CB, two AirNow stations (Horn Point and Hampton) near the south of CB, and one AirNow station (Albemarie_HS) to the west of Atmosphere 2022, 13, 641 9 of 20 CB. As expected, the hourly surface O 3 mixing ratio diurnal variation showed a lower O 3 mixing ratio at night (around 10 ppbv at 0600 UTC; sunrise was at 1030 UTC), and a higher O 3 mixing ratio in the daytime (around 30 ppbv at 2000 UTC; sunset was at 0000 UTC). During the night, the absence of sunlight and fewer O 3 precursor (NO X , etc.) emissions prevented the required O 3 -forming chemical reactions from being initiated. The correlation coefficients (Rs) and root mean square errors (RMSEs) were 0.77, 0.95, 0.93 and 6.38, 7.99, 9.31, 0.82, 0.65, and 0.91 for Milington, Fair Hill, Aldino, Horn Point, Hampton, and Albemarie_HS, respectively. In general, the WRF-Chem model captured the AirNow-measured O 3 mixing ratio diurnal variation trend pattern but overestimated the O 3 concentration by 5-10 ppbv, or 20-30%.

Overview of the O 3 Mixing Ratio Difference
The observational and modeled O 3 mixing ratio diurnal variation revealed the transition between the minimum and maximum of O 3 mixing ratios ( Figure 6). By separating the nighttime modeled O 3 mixing ratio difference from that of the daytime, the absolute and relative influences of CB can be analyzed explicitly. Daytime is defined as 1000-0000 UTC (or 0600-2000 local time), and nighttime is defined as 0000-0900 UTC (or 00:00-06:00 and 20:00-00:00 local time). The relative difference is calculated using the following equation.
On 3 June, at night, the model-simulated O 3 with water over CB was 28 ppbv, while, without water, the model-simulated O 3 mixing ratio was 18 ppbv (Figure 7a,b). In the day, the model-simulated O 3 with water over CB was 32 ppbv, while, without water, the model-simulated O 3 mixing ratio was 25 ppbv (Figure 7c,d). The averaged surface O 3 mixing ratio over CB increased by 10 ppbv (30%) at night (Figure 7c,d) and 5 ppbv (20%) during the day (Figure 7g,h). The shallow stable nocturnal boundary layer "trapped" O 3 and prevented it from venting to higher altitudes at night, while the deep mixed layer lofted O 3 to the free troposphere during the daytime. In the daytime, the O 3 mixing ratio difference was distributed semi-homogeneously over CB, with more O 3 over southern CB than northern CB (Figure 7g,h). The higher O 3 increase over southern CB was due to the larger area of the water body in southern CB. The O 3 mixing ratio difference increased slightly from east to west across CB, consistent with the direction of the prevailing winds. The wind difference over CB indicated wind increase due to less surface friction over the water surface than over land and contributions from bay breezes. In addition, the westward pooling of pollution created elevated O 3 that extended into the Baltimore-Washington D.C. urban corridor [26]. At night, the O 3 mixing ratio difference gradient gradually decreased and approached the minimum close to the eastern shore (Figure 7c,d). In the daytime, the O 3 mixing ratio difference gradient had a pattern analogous to nighttime but had a smaller magnitude (Figure 7g,h). The surface O 3 mixing ratio diurnal variations in both the baseline and sensitivity experiments and their differences above CB are shown in Figure 8. These show that the surface O 3 mixing ratio above CB from the baseline and sensitivity experiments had a similar diurnal variation trend. The surface O 3 mixing ratio at night (about 20-25 ppbv) was lower than in the daytime (about 25-30 ppbv). The larger surface O 3 mixing ratio difference mean over CB occurred at night (about 5 ppbv or 25%) rather than in the daytime (about 2 ppbv or 15%) due to the longer residence time of O 3 over CB at night. The wind speed over CB at night is lower than in the daytime, recorded as 8 m s −1 and 12 m s −1 , respectively. At night, the residence of this higher O 3 air over CB tended to be longer. In addition, O 3 accumulated at the southern end of CB (and the western shore of CB) where it was more impacted, and differences would add to the chance of exceedance of O 3 mixing ratio limits. In fact, Goldberg et al. (2014) indicated that over a 10-day period, ozone was exceeded four times over the water as observed on a moving ship compared to two times over land during the DISCOVER-AQ-2011 project (June-31 July 2011, [26]), while the surrounding land exceeded the NAAQS four (4) times, reinforcing the O 3 increase effect due to CB, especially in the downwind region.

Dynamical Influence on O 3 Mixing Ratios 4.2.1. Horizontal Dynamical Influence
We analyzed the dynamic mechanism that increased the simulated O 3 mixing ratio solely due to the change of the land surface type. First, the thermal properties of water and land (loam) are different. The volumetric heat capacity of land (loam) is around 2.0 × 10 6 J m −3 K −1 , much smaller than water, 4.18 × 10 6 J m −3 K −1 . This implies that the loam temperature increases/decreases about twice as much as water by absorbing/emitting the same amount of heat during the day/night. This fundamental physical property controls the model lower boundary surface temperature difference in our experiments. Due to the presence of a water body over CB, the temperature difference increased by 2.0-3.0 K at night, while the temperature difference decreased by 1.0 K in the day at the northern end of CB (Figure 9a,b). The nighttime temperature difference increase over CB was distributed homogeneously across northern and southern CB. However, the daytime temperature change showed a regional variation, the temperature over northern CB decreasing by about 1 K compared to a decrease of only 0.1 K or even an increase of 0.1 K over southern CB. The absolute value of temperature increase was larger than the temperature decrease. The obvious reason was that land (loam) emits longwave radiation and its temperature decreases quickly, whether it is night or day. The absence of the shortwave radiation at night led to the temperature decrease over land being relatively large. During the day, besides the longwave radiation emission, there was also shortwave radiation absorption.
Land increased temperature faster than water, while absorbing the same amount of heat. So, the water temperature was slightly lower during the day. The regional temperature variation also led to water vapor mixing ratio variability over CB. The resulting water vapor mixing ratio differences increased both at night (~2 g kg −1 ) and in the day (~0.5 g kg −1 ) (Figure 9c,d). Similarly, the vertical velocity difference also changed from night to day (Figure 9e,f). The vertical velocity changes were located at western CB and also mostly at northern CB. There was a velocity increase of about 0.005 m s −1 at the western shore of northern CB at night and a decreased of 0.002 m s −1 at the same location during the day, which corresponded to the temperature difference shown in Figure 9a,b. The positive water surface temperature change led to air ascent at night, while the negative water surface temperature change caused air to descend. The experiment showed that the water body modified the meteorology and dynamics of CB in a major way, leading to mixing and advection (via the prevailing winds). The higher surface temperature led to instability near the surface and more water vapor flux through the relatively higher vertical velocity, hence the flux that resulted. Figure 9. WRF-Chem-simulated absolute differences (water minus nowater) in temperature (a,b), water vapor mixing ratio (c,d), and vertical velocity (e,f) between simulations of baseline (water) and sensitivity (nowater) experiments over CB overlaid with corresponding surface wind differences. The first row is at night; the second row is in the day. Daytime is defined as 1000-0000 UTC (or 06:00-20:00 local time); nighttime is defined as 0000-0900 UTC (or 00:00-06:00 and 20:00-00:00 local time).
Pressure differences were observed from the baseline and sensitivity experiments. At night, the effect of CB increased the surface pressure over the western shore of northern CB, leading to the pressure gradient force pointing from west to east (Figure 10a,c). This is consistent with Goldberg et al. (2014) [26]. The pressure gradient contributed to the land breeze from the west, offshore of CB. The resulting increased land breeze blew in opposition to the northeasterly prevailing wind, minimizing the wind difference over CB due to the water body (Figure 10a,c). During the daytime, CB decreased pressure. The pressure gradient favored the prevailing wind and worked to reinforce the flow (Figure 10b,d). The differences in surface temperature, moisture, and vertical velocity have important influences on the stability and dynamics of the boundary layer. To investigate the stability consequences of these meteorological quantities, stability parameters were investigated. First, the gradient Richardson number (Ri) was defined as: The gradient Richardson number is a dynamical stability measure to determine whether turbulence exists. g = 9.8 m s −2 is the gravitational acceleration. U and V are horizontal velocities. θ v is the virtual potential temperature and can be calculated as: r is the water vapor mixing ratio, and r L is the liquid water mixing ratio in the atmosphere. For the unsaturated atmosphere, r L = 0, and ( ∂U ∂Z ) 2 + ∂V ∂Z 2 is the wind shear squared.
The root of the numerator in Equation (1) is the Brunt-Väisälä frequency (N), a measure of buoyancy, calculated using the following equation: If the Brunt-Väisälä frequency is positive, the atmosphere is stable; if the Brunt-Väisälä frequency is negative, the atmosphere becomes unstable. Figure 11 shows buoyancy, Brunt-Väisälä (buoyancy) frequency, wind shear squared, and gradient Richardson number at the surface, with a relative decrease in magnitude and no regional preference at night, which led to more instability due to the existence of CB. However, during the day, positive values for all parameters were found located over northern CB, while negative values were found at southern CB.

Vertical Dynamical Influence
Apart from the O3 mixing ratio horizontal distribution, its vertical distribution also helps to understand how CB influences the O3 three-dimensional distribution. The cross sections should pass through the east-west, north-south, and center of CB. We selected six cross sections to visualize the vertical pictures around CB. The vertical cross sections of O3 mixing ratio and water vapor mixing ratio were selected to characterize the role of CB in impacting their vertical distributions. Figure 12 shows the O3 mixing ratio vertical distribution difference throughout northern, central, and southern CB to distinguish the possible dynamic variations due to At night, the buoyancy, Brunt-Väisälä (buoyancy) frequency, wind shear squared, and gradient Richardson number decreased due to the presence of CB (Figure 11a,c,e,g). A relatively small water body surrounded by a large landmass over northern CB led to a significant thermally induced perturbation in contrast to southern CB. Relatively small surface wind differences were due to the water body and contributed to the atmospheric stability through wind shear. There was no regional preference of surface horizontal wind difference and thus its contribution to the instability was uniform over the entirety of CB (Figure 11a,c,e,g). In addition to the surface temperature increase, the vertical velocity difference was positive over northern and middle CB, which contributed to the unstable condition (Figure 10e). The critical Richardson number is 0.25. If the gradient Richardson number is less than 0.25, the atmosphere is dynamically unstable and turbulent. The gradient Richardson number difference was negative (around −1-−2) at northern CB and caused significant atmospheric instability. In the day, the surface wind differences at northern CB were smaller than at nighttime. The differences in buoyancy, Brunt-Väisälä (buoyancy) frequency, wind shear squared, and gradient Richardson number remained negative at southern CB and shifted to positive values at northern CB (Figure 11b,d,f,h). Contrary to the night, CB increased the gradient Richardson number in the day at northern CB, leading to the stable condition.

Vertical Dynamical Influence
Apart from the O 3 mixing ratio horizontal distribution, its vertical distribution also helps to understand how CB influences the O 3 three-dimensional distribution. The cross sections should pass through the east-west, north-south, and center of CB. We selected six cross sections to visualize the vertical pictures around CB. The vertical cross sections of O 3 mixing ratio and water vapor mixing ratio were selected to characterize the role of CB in impacting their vertical distributions. Figure 12 shows the O 3 mixing ratio vertical distribution difference throughout northern, central, and southern CB to distinguish the possible dynamic variations due to different parts over CB. The PBLHs are shown for both the baseline and sensitivity experiments. At night, O 3 was constrained within the shallow nocturnal boundary layer and the PBLH was less than 0.4 km, which led to the O 3 increase at night as compared to the daytime (Figure 12a-c). The increased O 3 mixing ratio difference was found near the surface to increasingly intrude over the western shore via transport by the nighttime breeze and the northeasterly prevailing wind. The water body of CB increased the PBLH by about 0.1 km over the west shore of northern CB. A subsidence of about 0.1 km in PBLH over eastern CB was also evident due, presumably, to the resulting subsidence (Figure 12a). Similar dynamics were observed at the cross sections over middle and southern CB (Figure 12b,c). After sunrise, the resulting instability and convection increased the PBLH from 0.4 to 1.0 km. PBLHs were tilted from west to east of CB. The west is close to land terrain, which meant a higher temperature during the day, while the east is near the Atlantic Ocean which had a lower temperature during the day. All these factors diluted the O 3 mixture and transported O 3 to the free troposphere, leading to surface O 3 mixing ratio difference decrease.
Another interesting point is the strength of the wind circulation in the vertical cross sections shown in Figure 12. The wind vectors for baseline and sensitivity experiments are not shown here due to the large number of figures. As the bay breeze moved from north to center to south, it became well defined and stronger (shown as circulations in Figure 12e,f). At the northern CB transect, surrounded by land and far from the Atlantic Ocean, the simulated PBLH increased dramatically when the water body was replaced by land (loam), implying that the land surface contribution to the PBLH was large. This indicated an ability to absorb large amounts of heat (shortwave radiation) with relatively little temperature change and resulted in subsiding air over CB. The overall result was that the maximum O 3 mixing ratio difference was localized over western CB due to the shallow PBL with a relatively weak bay breeze to ventilate the O 3 buildup. The change in PBLH was smaller in the central and southern transects where a much stronger breeze and shallower PBLH were observed. factors diluted the O3 mixture and transported O3 to the free troposphere, leading to surface O3 mixing ratio difference decrease. Figure 12. WRF-Chem-simulated O3 mixing ratio vertical cross section (from west to east) differences (water minus nowater, night: (a-c); day: (d-f)). The overlaid streamline is the corresponding surface wind difference. The green curve is the PBLH of water; the pink curve is the PBLH of nowater; CB is represented by the yellow bar. The upright insert map shows the location of the vertical cross section (blue line). Note that vertical wind is magnified by 50 times for the illustration. Daytime is defined as 1000-0000 UTC (or 06:00-20:00 local time); nighttime is defined as 0000-0900 UTC (or 00:00-06:00 and 20:00-00:00 local time).
Another interesting point is the strength of the wind circulation in the vertical cross sections shown in Figure 12. The wind vectors for baseline and sensitivity experiments are not shown here due to the large number of figures. As the bay breeze moved from north to center to south, it became well defined and stronger (shown as circulations in Figure  12e,f). At the northern CB transect, surrounded by land and far from the Atlantic Ocean, the simulated PBLH increased dramatically when the water body was replaced by land (loam), implying that the land surface contribution to the PBLH was large. This indicated an ability to absorb large amounts of heat (shortwave radiation) with relatively little temperature change and resulted in subsiding air over CB. The overall result was that the maximum O3 mixing ratio difference was localized over western CB due to the shallow PBL with a relatively weak bay breeze to ventilate the O3 buildup. The change in PBLH was smaller in the central and southern transects where a much stronger breeze and shallower PBLH were observed. Figure 13 shows O3 mixing ratio north-south transects at the west, over, and east of CB, which confirms the significant influence of the water surface and associated bay dynamics on the O3 mixing ratio and its distribution over and west of CB (Figure 13a-d), while no influence is visible at the east of CB (Figure 13e,f). The northeastern prevailing wind blew O3 to the west after its generation from CB. Figure 13 also shows the shallow and flat nocturnal PBLH (~0.2 km) at night as compared to the tilted and large PBLH in the day. We did not find the significant PBLH difference between water and nowater simulations at nighttime, whereas a totally different scenario unfolded for daytime PBLH over CB ( Figure  13d), while the water-influenced PBLH was similar to the sensitivity simulations both west and east of CB, decreasing sharply from north to south, perhaps due to the water body area difference between northern and southern CB (Figure 13b,f). The PBLH over CB from the baseline experiment was much lower than the sensitivity experiment, approximately 0.8 km Figure 12. WRF-Chem-simulated O 3 mixing ratio vertical cross section (from west to east) differences (water minus nowater, night: (a-c); day: (d-f)). The overlaid streamline is the corresponding surface wind difference. The green curve is the PBLH of water; the pink curve is the PBLH of nowater; CB is represented by the yellow bar. The upright insert map shows the location of the vertical cross section (blue line). Note that vertical wind is magnified by 50 times for the illustration. Daytime is defined as 1000-0000 UTC (or 06:00-20:00 local time); nighttime is defined as 0000-0900 UTC (or 00:00-06:00 and 20:00-00:00 local time). Figure 13 shows O 3 mixing ratio north-south transects at the west, over, and east of CB, which confirms the significant influence of the water surface and associated bay dynamics on the O 3 mixing ratio and its distribution over and west of CB (Figure 13a-d), while no influence is visible at the east of CB (Figure 13e,f). The northeastern prevailing wind blew O 3 to the west after its generation from CB. Figure 13 also shows the shallow and flat nocturnal PBLH (~0.2 km) at night as compared to the tilted and large PBLH in the day. We did not find the significant PBLH difference between water and nowater simulations at nighttime, whereas a totally different scenario unfolded for daytime PBLH over CB (Figure 13d), while the water-influenced PBLH was similar to the sensitivity simulations both west and east of CB, decreasing sharply from north to south, perhaps due to the water body area difference between northern and southern CB (Figure 13b,f). The PBLH over CB from the baseline experiment was much lower than the sensitivity experiment, approximately 0.8 km lower (Figure 13d). This was explained, as previously discussed, by the significant impact of the much smaller water body area over northern CB.
Water vapor can influence atmospheric stability as well. Equation (2) shows that the calculation of the gradient Richardson number uses the virtual potential temperature, which is influenced by water vapor (Equation (3)). If the water vapor increases near the surface, the virtual potential temperature will increase, which leads to a decrease in the virtual potential temperature vertical gradient. Then, the gradient Richardson number will decrease and the atmosphere will become unstable. The water body in CB influenced the water vapor mixing ratio vertical distribution (Figure 14). At night, the presence of CB increased the water vapor mixing ratio from north to south CB (Figure 14a,c,e). The water vapor increased much less over northern CB than over southern CB (Figure 14a,e), while, during the day, the water vapor mixing ratio increased less than at nighttime (Figure 14d,f) or even decreased over northern CB (Figure 14b). The possible reason is that the temperature decreased (Figure 9b). Overall, the atmosphere became unstable due to the increases in water vapor over CB. lower (Figure 13d). This was explained, as previously discussed, by the significant impact of the much smaller water body area over northern CB. Figure 13. WRF-Chem-simulated O3 mixing ratio vertical cross section (from north to south) differences (water minus nowater, night: (a,c,e); day: (b,d,f)). The overlaid streamline is the corresponding surface wind difference. The green curve is the PBLH of water; the pink curve is the PBLH of nowater. CB is denoted with the yellow bar. Note that vertical wind is magnified by 50 times for illustration purposes. Daytime is defined as 1000-0000 UTC (or 06:00-20:00 local time); nighttime is defined as 0000-0900 UTC (or 00:00-06:00 and 20:00-00:00 local time).
Water vapor can influence atmospheric stability as well. Equation (2) shows that the calculation of the gradient Richardson number uses the virtual potential temperature, which is influenced by water vapor (Equation (3)). If the water vapor increases near the surface, the virtual potential temperature will increase, which leads to a decrease in the virtual potential temperature vertical gradient. Then, the gradient Richardson number will decrease and the atmosphere will become unstable. The water body in CB influenced the water vapor mixing ratio vertical distribution (Figure 14). At night, the presence of CB increased the water vapor mixing ratio from north to south CB (Figure 14a,c,e). The water vapor increased much less over northern CB than over southern CB (Figure 14a,e), while, during the day, the water vapor mixing ratio increased less than at nighttime (Figure 14d,f) or even decreased over northern CB (Figure 14b). The possible reason is that the temperature decreased (Figure 9b). Overall, the atmosphere became unstable due to the increases in water vapor over CB. Figure 13. WRF-Chem-simulated O 3 mixing ratio vertical cross section (from north to south) differences (water minus nowater, night: (a,c,e); day: (b,d,f)). The overlaid streamline is the corresponding surface wind difference. The green curve is the PBLH of water; the pink curve is the PBLH of nowater. CB is denoted with the yellow bar. Note that vertical wind is magnified by 50 times for illustration purposes. Daytime is defined as 1000-0000 UTC (or 06:00-20:00 local time); nighttime is defined as 0000-0900 UTC (or 00:00-06:00 and 20:00-00:00 local time).
Atmosphere 2022, 13, x FOR PEER REVIEW 18 of 21 Figure 14. WRF-Chem-simulated water vapor mixing ratio vertical cross section (from west to east) differences (water minus nowater, night: (a-c); day: (d-f)). The overlaid streamline is the corresponding surface wind difference. The green curve is the PBLH of water; the pink curve is the PBLH of nowater; CB is represented by the yellow bar. The upright insert map shows the location of the vertical cross section (blue line). Note that vertical wind is magnified by 50 times for the illustration. Daytime is defined as 1000-0000 UTC (or 06:00-20:00 local time); nighttime is defined as 0000-0900 UTC (or 00:00-06:00 and 20:00-00:00 local time).

Discussion and Conclusions
This paper has investigated the influence of CB on local O3 pollution by employing the WRF-Chem model to simulate the O3 concentration on 3 June 2015. To ensure model accuracy, we first validated the model's performance in a simulation of O3 spatial distribution and time evolution using AirNow surface O3 observations. Through the comparison between baseline Figure 14. WRF-Chem-simulated water vapor mixing ratio vertical cross section (from west to east) differences (water minus nowater, night: (a-c); day: (d-f)). The overlaid streamline is the corresponding surface wind difference. The green curve is the PBLH of water; the pink curve is the PBLH of nowater; CB is represented by the yellow bar. The upright insert map shows the location of the vertical cross section (blue line). Note that vertical wind is magnified by 50 times for the illustration. Daytime is defined as 1000-0000 UTC (or 06:00-20:00 local time); nighttime is defined as 0000-0900 UTC (or 00:00-06:00 and 20:00-00:00 local time).

Discussion and Conclusions
This paper has investigated the influence of CB on local O 3 pollution by employing the WRF-Chem model to simulate the O 3 concentration on 3 June 2015. To ensure model accuracy, we first validated the model's performance in a simulation of O 3 spatial distribution and time evolution using AirNow surface O 3 observations. Through the comparison between baseline and sensitivity experiments over CB, by altering the surface type of CB from water to land (loam), dynamical differences due to CB were analyzed thoroughly from both horizontal and vertical directions. The findings showed that CB elevated the O 3 mixing ratio both at nighttime and in the daytime, but for discrepant reasons.
The model-simulated O 3 mixing ratio overall agreed with AirNow O 3 observations, even though it overestimated ratios near mountains and underestimated them over New Jersey and Delaware. The higher O 3 mixing ratio over CB was decidedly confirmed from the perspective of the model simulation, which was consistent with the previous studies over CB [24,26,27]. A shallower boundary layer decreased the surface O 3 vertical dilution, especially for northern CB, which led to relatively high O 3 concentration close to the surface. A strong bay breeze caused O 3 to accumulate in the downwind region near CB. Taking the model uncertainty into account, the surface O 3 mixing ratio increased by 10% at night and 5% in the day. The O 3 mixing ratio difference gradient was negative over the eastern shore and positive over the western shore, which indicated that the O 3 mixing ratio difference increased from the eastern shore to the western shore due to the northeastern prevailing winds. With the difference in O 3 mixing ratio diurnal variation between the baseline and sensitivity experiments, water increased O 3 by 5 ppbv (25%) at night and 2 ppbv (15%) during the day. The water body increased temperature by 2-3 K at night and decreased it by 1.0 K during the day, increased the water vapor mixing ratio by 2 g kg −1 at night and 0.5 K in the day, and increased vertical velocity by 0.005 m s −1 at night and decreased it by 0.002 m s −1 .
The vertical distributions of the O 3 mixing ratio and the water vapor mixing ratio showed the vertical influence of CB. The west-east cross sections located at northern, central, and southern CB showed that O 3 was confined within PBL, which reached 0.2 km at night and extended to 0.8-1.0 km in the day (tilted PBL). PBLH is lowest over CB. Over land, PBLH is lower close to the upwind cold water, i.e., the eastern shore, and downwind and close to the Atlantic Ocean. PBLH is higher on the western shore of CB, where an urban heat island is present.
Apart from field campaigns, this work sheds light on how to employ models to study the influence of CB on O 3 mixing ratios thoroughly, since models provide us with the flexibility to eliminate CB. This helps us to quantify the O 3 concentration differences between them, with and without CB, even though model error attributes some uncertainties to these estimations. A merit of this modeling approach is that it provides us with a novel method to study the influence of CB on O 3 in three dimensions, which could not be achieved with field campaigns. As models are further developed, the model errors will become smaller and analogous work in the future will be more reliable.  Data Availability Statement: MODIS products are publicly available at: https://atmosphereimager.gsfc.nasa.gov/products.html, accessed on 10 May 2020. The EPA AirNow measurements are publicly available at: https://www.airnowtech.org/index.cfm, accessed on 10 June 2020. The WRF-Chem model is publicly available at: https://www2.mmm.ucar.edu/wrf/users/, accessed on 2 February 2020.