WRF Cold-Pool Sensitivities Relevant to Winter-Ozone Modelling in the Uinta Basin

1. Introduction

Surface wintertime ozone in the Uinta Basin develops within a snow-covered, persistently stable boundary layer in which weak mixing and transport can confine precursors near the surface [1,2,3]. For a coupled weather–chemistry calculation, the relevant meteorological information therefore extends beyond 2 m temperature error to the simulated depth, stability, heat storage, and transport of the lower atmosphere. Point observations remain indispensable, but a horizontally separated surface network does not by itself validate inversion depth or the three-dimensional mixing volume supplied to chemistry.

Meteorological, emissions, photolysis, deposition, and chemical-mechanism errors can compensate for one another in coupled simulations. For example, holding meteorology fixed while changing emissions can alter ozone agreement without identifying a unique inventory error [1]. Cumulative changes in emissions, leaf area, deposition, and meteorology likewise complicate attribution, and improved modelled temperature profiles need not translate directly into improved ozone profiles [4]. Credible meteorology is thus a prerequisite for chemical attribution rather than proof that meteorology is the dominant or only uncertainty.

Previous WRF studies also show why the experimental design must remain explicit. Bundled snow, microphysics, and grid choices can produce large case-specific responses, while evaluation and chemistry grids need not sample the same structure [2]. Four-dimensional data-assimilation responses can depend on the nudged variables, input data, and scheme-dependent boundary-layer diagnosis [5]. The present study therefore treats controlled contrasts, failed tests, and within-run diagnostics as distinct kinds of evidence.

Following preliminary experiments and review of the archived simulation family, we ask:

1.

Does the lower-layer response to the supplied GFS or NAM driving analysis recur under both one-way and two-way nest feedback?

2.

How large are the controlled nest-feedback, boundary-layer and surface-layer, and slope-aware-radiation responses relative to the driving-product response on common diagnostics?

3.

What reproducible atmospheric response appears across controlled static-terrain-source brackets on the same WRF grid, and which interpretations remain unsupported?

4.

What do surface, common-height, and footprint heat-deficit diagnostics reveal together that any single station metric can miss?

Terrain-source contrasts require area-weighted geometry diagnostics because a station-only summary can miss a rim- or slope-confined response. The archived simulations use the compact X0–X10 identifiers, which link configurations in the text and tables to the run inventory. The present paper evaluates the meteorological state that a later chemistry calculation would receive; it contains no chemistry calculation and makes no ozone prediction.

2. Materials and Methods

2.1. Case, Observations, and Evaluation Domains

The case is 2 February 2013 within a broader January–February period of snow cover, cold-pool conditions, and elevated ozone in the Uinta Basin [1,2,5]. The WRF integrations cover 1200–1800 Coordinated Universal Time (UTC), corresponding to 0500–1100 Mountain Standard Time (MST; UTC$-7$). The retained analysis window covers 1400–1800 UTC, or 0700–1100 MST, and excludes both the 1200 UTC initialised output and the 1300 UTC first-hour forecast. This leaves five hourly samples for every common-window contrast. Surface observations were obtained through the Synoptic Data PBC Mesonet API for 2 m temperature at 11 sites, 10 m wind speed at 10 sites, and surface pressure at 9 sites.1 Subhourly reports were mean-averaged within each UTC clock hour. The stations span basin-floor, sidewall, and bench locations but are neither spatially uniform nor area weighted; as the focus herein is not verification against observations, we leave further evaluation to future studies. The 1 km d02 domain was selected for the common observation comparisons because seven stations fall outside the usable d03 interior, leaving four usable d03 station columns. Figure 1 combines the reporting-station geometry with the full nested-domain layout. Use of d02 does not establish that d02 and d03 produce the same response, and a common-support transfer test remains future work. Four stations spanning the basin floor and benches provide surface context from 31 January 0000 UTC to 4 February 0000 UTC. Observed and modelled surface potential temperature for a two-site structural index was calculated as $\theta =T{(1000/p)}^{0.2854}$, with T in kelvin and p in hPa. The index is Mountain Home on the north-west bench minus Pariette Draw on the basin floor, with both locations labelled in Figure 1. Their horizontal separation requires treating the index as a surface contrast rather than free-air inversion strength or cold-pool depth.

Episode-scale ozone context came from the US Environmental Protection Agency (EPA) Air Quality System (AQS) AirData daily archive for 1 January–28 February 2013. We retained the EPA-reported “Ozone 8-hour 2015” daily-summary row at Ouray, Redwash, Roosevelt, Vernal, and Little Mountain. All five sites contained one daily row on each of the 59 local calendar dates; within-day observation coverage ranged from 47% to 100% and is marked explicitly in the figure. The satellite panel in Figure 2 uses the National Aeronautics and Space Administration (NASA) Global Imagery Browse Services (GIBS) Terra Moderate Resolution Imaging Spectroradiometer (MODIS) corrected-reflectance 7-2-1 layer for 2 February 2013.2

2.2. WRF Configuration and Archived Simulation Family

The controlled six-hour simulations used WRF version 4.8.0, initialised at 1200 UTC with either Global Forecast System (GFS) or North American Model (NAM) analyses and supplied with source-specific lateral boundary conditions through 1800 UTC on 2 February 2013. These are retrospective simulations driven by analyses rather than six-hour operational forecasts. The design does not separate the initialised state from the subsequent lateral-boundary contribution. We used three nested computational domains, d01–d03, with horizontal grid spacings of 3 km, 1 km, and 333.333 m. The controlled family used 75 eta levels concentrated near the surface, with approximately 20 m separation among the lowest mass levels, and a 50 hPa model top. The Advanced Research WRF Model Version 4 is cited following Skamarock et al. (2019, NCAR/TN-556+STR).3 The reference physics suite used Thompson microphysics, Rapid Radiative Transfer Model for General Circulation Models (RRTMG) longwave and shortwave radiation, the revised fifth-generation Pennsylvania State University–National Center for Atmospheric Research Mesoscale Model (MM5) similarity surface layer, the Noah land-surface model, and the Yonsei University (YSU) planetary boundary layer (PBL) scheme.4 Grid and observation four-dimensional data assimilation were disabled. Departures from that reference configuration are listed by simulation in Table 1.

Static terrain entered through the WRF Preprocessing System and was held fixed within each controlled pair. The custom fine source used 3-arc-second terrain, the default source was 30-arc-second class, and the deliberately coarse source used the 5-arc-minute GMTED2010 product. At approximately ${40}^{\circ }$N their nominal north–south spacings are about 90 m, 900 m, and 9 km, respectively, and these static-source spacings are distinct from WRF horizontal grid spacing. Land-surface initial fields, including snow, entered through the source-specific meteorological preprocessing, so the GFS–NAM treatment does not isolate a snow-field contribution. Our experiment configurations split into four scientific branches: driving-data product, boundary-layer and radiation treatments, nest feedback, and static-terrain source. The X7 and X8 configurations differ from the controlled family in several dimensions at once and enter only as design context (Section 2.4) — again, interpretation must be careful when changes in output may result from confounding effects.

2.3. Controlled Contrasts and Estimands

We contrast simulations, identifying the dominant lever on the simulation’s sensitivity, and its magnitude (Table 2). The table is ordered by scientific role. The four default-terrain simulations X3, X5, X6, and X10 form a fundamental two-by-two matrix that buttresses our analyses. Native GFS and NAM data provides 27 and 40 metgrid levels, respectively.

(The remaining controlled contrasts enter on their own design and evidence: the X2-minus-X1 scheme treatment, X1-minus-X0 slope-radiation treatment, X9-minus-X0 fine-terrain feedback treatment, and the pairwise X3/X0/X4 terrain-source brackets.)

2.4. Design Exclusions: Resolution and Legacy Context

X7, X0, and X8 do not form a controlled resolution sequence because their integrations differ in forecast age, domain geometry, feedback, vertical levels, terrain placement, radiation settings, and output cadence. They are retained only to document the exclusion; no X8 diagnostic is reported in this paper.

2.5. Surface and Common-AGL Diagnostics

Modelled 2 m dry-bulb temperature, 10 m wind, and surface pressure were sampled at the nearest d02 cell to each station. None of the selected cells lay on steep terrain or an obviously unrepresentative model location. Across the eleven fixed columns, model terrain offset from reported station elevation spanned -32 m to 14 m. While this justifies the nearest-neighbour approach, the pairing of WRF output to observation value does not remove (wind, solar) exposure, slope angle, land-cover use, or network-coverage deficiency. Mass-level height above ground level (AGL) was calculated from perturbation plus base geopotential relative to local WRF terrain height. Grid-relative winds were destaggered and rotated to earth-relative components. Model-level potential temperature $\theta$ and horizontal wind components were interpolated linearly at 25 m increments from 50 to 500 m AGL at each station column. Fixed-layer stability was defined as

$$S_{50:500}=\theta \left(500\mathrm{m}\right)-\theta \left(50\mathrm{m}\right)$$

(1)

The transport-wind proxy was defined as

$$U_{50:500}={\left({\overline{u}}^2+{\overline{v}}^2\right)}^{1/2}.$$

(2)

where the bars denote height means from 50 to 500 m AGL. This wind-derived quantity is a fixed-layer transport proxy and does not sample flow through a basin boundary, pollutant residence time, cold-pool volume, or a chemistry-model ventilation coefficient. Native WRF PBL height was retained only as a scheme-dependent companion diagnostic, consistent with the diagnostic cautions discussed by Tran et al. [5].

2.6. Pairing, Aggregation, and Short-Sample Sensitivity

Run differences were paired by station and valid hour before spatial aggregation. Each fixed station field was reduced to a single unweighted spatial median per hour; the five-hour value is thus a mean of those five medians. The stations therefore define one sampled spatial field and are not treated as independent replicates. A descriptive circular moving-block calculation drew 20,000 resamples of the five hourly medians with two-hour blocks and reports the 2.5th and 97.5th percentiles as a sensitivity interval. Because this is one deterministic case, interpretation emphasises practical effect size rather than population inference. The block ranges are descriptive dependence checks rather than confidence intervals.

2.7. Footprint-Integrated Heat Deficit

We calculated a crest-referenced atmospheric heat-deficit proxy following the basin-volume formulation of Whiteman et al. (1999) and the closest column implementation of Adler et al. (2023).5 Let $p_s$ and $p_c$ denote the pressures at the surface and fixed crest reference, respectively, with $p_c<p_s$. The present area-normalised atmospheric heat-deficit proxy was defined as $H=(c_p/g){\int }_{p_c}^{p_s}max\{{\theta }_{\mathrm{crest}}-\theta \left(p\right),0\}dp$, where ${\theta }_{\mathrm{crest}}$ is the potential temperature interpolated to the crest height in each modelled column. H is reported as the cell-area-weighted, footprint-mean, vertically integrated heat deficit in MJ m⁻². The primary footprint is the convex hull of eleven evaluation stations and eleven registered basin waypoints intersected with model terrain below the 2200 m crest reference. It contains 4,965 d02 cells and covers 5,277 km². The seven predeclared definitions were the combined station-and-waypoint hull at crest references of 2100, 2200, and 2300 m; the 2200 m combined hull dilated or eroded by one cell; and separate station-only and waypoint-only hulls at 2200 m. Each feedback-matched GFS–NAM pair used a common terrain mask. The weaker-state criterion required the GFS-minus-NAM difference to remain negative at every retained hour under all seven definitions. The faster-loss criterion required both a negative GFS change from 1400 to 1800 UTC and a more negative change than in NAM under all seven definitions. Both criteria were fixed before output inspection, and 1300 UTC was retained only as trajectory context. The footprint is an evaluation region rather than a hydrologic basin boundary, and the proxy is neither a measured basin-volume heat deficit nor an observed cold-pool quantity.

3. Results

3.1. Observed Episode Context

The MODIS browse image shows extensive snow and ice across the Uinta Basin and surrounding terrain on 2 February (Figure 2a). The AQS series places the simulation date within a broader period of elevated daily ozone maxima (Figure 2b). On 2 February, the reported daily maximum eight-hour values were 56 ppb at Little Mountain, 65–70 ppb at Vernal and Roosevelt, 81 ppb at Redwash, and 95 ppb at Ouray. These observations provide episode context and do not establish meteorology–ozone attribution.

3.2. Observed Surface-Temperature Context

The basin-floor station repeatedly separated from the higher sites between 31 January and 3 February, while the magnitude and ordering of station temperatures evolved through the diurnal cycle (Figure 3). The recurring floor-to-bench separation is consistent with a cold-pool surface pattern. Because the stations are horizontally separated, however, the time series do not determine the cold-pool top, depth, volume, or vertical temperature structure.

3.3. Common-Pipeline Contrasts at a Glance

Figure 4 compares every controlled contrast evaluated through the common station and common-AGL pipeline. The two driving-product contrasts are farthest from zero for temperature, fixed-layer stability, and transport wind. The PBL/surface-layer treatment is smaller, while slope radiation, nest feedback, and the driving-product-by-feedback interaction cluster near zero on those same scales. Pressure and wind-error metrics follow different orderings, so the panel is a response hierarchy rather than an overall configuration ranking.

3.4. Driving-Product Response in the Common Lower Layer

The two-by-two matrix in Figure 5 separates the GFS-versus-NAM driving-product response under one-way and two-way feedback. Relative to GFS, NAM was 3.63 K more stable with two-way feedback and 3.65 K more stable with one-way feedback through the 50–500 m layer (Table A1). The five-hour means of hourly spatial-median stability were 7.50 and 7.38 K for the two GFS runs, compared with 10.82 K for both NAM members. The hourly GFS-minus-NAM differences ranged from $-5.58$ to $-1.13$ K with two-way feedback and from $-5.36$ to $-0.90$ K with one-way feedback. All ten product-by-hour stability contrasts therefore retained the same sign. GFS also increased fixed-layer transport wind by 0.79 m s⁻¹ with two-way feedback and 0.99 m s⁻¹ with one-way feedback. The corresponding five-hour means of hourly spatial-median transport wind were approximately 2.0 m s⁻¹ for both GFS members and approximately 1.1 m s⁻¹ for both NAM members. The hourly differences ranged from 0.37 to 1.18 m s⁻¹ with two-way feedback and from 0.61 to 1.50 m s⁻¹ with one-way feedback, with the same sign in every hour. Native PBL-height differences changed sign within both feedback modes and did not reproduce the stability and wind separation. This divergence is consistent with the scheme dependence of native PBL-height diagnosis discussed by Tran et al. [5].

3.5. Surface and Spatial Corroboration

The surface-temperature contrasts had the same product ordering as the common-layer diagnostics. GFS-minus-NAM 2 m temperature-bias contrasts were +1.79 °C with two-way feedback and +1.76 °C with one-way feedback (Table A2). Surface-pressure absolute-error contrasts changed in the opposite direction, at $-13.6$ and $-10.9$ Pa, while hourly wind absolute-error differences changed sign. No single surface score therefore provides a general product ranking. Figure 6 gives a separate two-site surface-structure comparison for the observations and two-way anchor simulations.

The five-hour mean surface potential-temperature contrast between Mountain Home and Pariette Draw was 16.94 K in the observations, 17.03 K in NAM two-way, and 9.95 K in GFS two-way. The close NAM value corroborates only this horizontally separated surface index; it does not validate modelled inversion depth, vertical structure, or overall surface skill. Across the registered footprint, the mean GFS-minus-NAM near-surface response was positive and the vertical response remained concentrated below 500 m AGL under both feedback modes (Table A2). The predeclared spatial-coherence test nevertheless failed because only 64–67% of footprint cells had the expected sign at 1400–1500 UTC, below the required two-thirds share at every retained hour. The result is therefore shallow and mean-coherent, but not spatially uniform throughout the retained window. Later hours were more coherent, but changing an every-hour rule after inspecting those maps would be post hoc. The maps remain useful descriptive diagnostics without converting the failed test into a passed one.

3.6. Footprint Heat Deficit: A Weaker State Without Faster Loss

The GFS runs retained less footprint-mean, vertically integrated heat deficit at every retained hour under all seven definitions (Section 2.7; Table A3). Mean GFS-minus-NAM differences were $-1.23$ MJ m⁻² with two-way feedback and $-1.25$ MJ m⁻² with one-way feedback. Primary-definition five-hour means were approximately 3.0 MJ m⁻² for both GFS runs and approximately 4.2 MJ m⁻² for both NAM runs. Every run gained heat deficit from 1400 to 1800 UTC, so the weaker GFS state did not result from demonstrated faster heat-deficit loss during the retained window. Figure 7 shows the four increasing trajectories and the persistent separation between the two GFS and two NAM members.

3.7. Process Treatments on the Fine-Terrain Branch

The MYJ/Eta versus YSU/revised-MM5 scheme pair changes only the PBL and surface-layer schemes while holding the fine terrain and slope-aware-radiation configuration fixed. MYJ/Eta increased 2 m-temperature bias by 0.37 °C, reduced 10 m-wind absolute error by 0.145 m s⁻¹, and increased fixed-layer transport wind by 0.098 m s⁻¹, with the same sign in all five hours (Supplementary Table S1). These changes are smaller than the driving-product response on the shared lower-layer metrics. The opposing temperature- and wind-error changes do not support an overall scheme-skill ranking.

The slope-radiation addition produced negligible response at the basin-floor station columns. A spatially organised skin-temperature response appeared on the sidewalls, but the archived flat-plane shortwave field and missing surface heat fluxes could not test the mechanism directly. The diagnostic is retained in Supplementary Figure S1 and does not show that slope-aware radiation improved or degraded the cold pool.

3.8. Nest-Feedback Robustness Across Two Terrain Regimes

Within-product feedback effects were small for both GFS and NAM on the default-terrain branch. The fine-terrain NAM pair reproduced that character: two-way minus one-way feedback increased stability by 0.04 K and reduced transport wind by 0.01 m s⁻¹. The descriptive block-sensitivity intervals defined in Section 2.6 spanned zero for both metrics. The consistently signed 10 m-wind and surface-pressure changes were also small and do not support a feedback recommendation beyond this five-hour case.

3.9. Terrain-Source Brackets

The terrain-source branch holds NAM forcing, one-way feedback, WRF grid spacing, and the YSU/MM5 physics configuration fixed while changing only the static-terrain source. The three brackets are default minus fine, default minus coarse, and fine minus coarse. The default-minus-fine response was small on d02 and failed its predeclared terrain-response test. The wider coarse brackets produced larger spatial responses, but the summary changes with domain support and diagnostic (Supplementary Table S2). The registered terrain summaries therefore describe a reproducible model response without establishing a preferred terrain source, drainage pathway, upstream reservoir, or basin cold-air-volume change.

4. Discussion

4.1. The Main Result Is a Difference in the Supplied Lower-Layer State

The present study does not seek a universal product or configuration ranking. Its principal result is a repeated difference in the lower-layer state supplied by the two driving products under both feedback modes.

• The GFS-driven runs supplied a warmer, less stable, and faster-moving 50–500 m layer than feedback-matched NAM members.
• The heat-deficit calculation reproduced the same ordering under every registered footprint definition.
• Surface pressure and wind error followed different orderings, preventing the response hierarchy from becoming a general skill ranking.

The GFS state is more consistent with weaker confinement than the analogous NAM state, but the transport-wind proxy and heat-deficit values are not substitutes for measured mixing volume or basin export. This distinction matters because errors in confinement, emissions, snow and cloud processes, photolysis, deposition, and chemical mechanism can compensate for one another in coupled simulations [1,2,4]. A practical future ensemble should therefore vary the treatments that show the largest phenomenon-relevant response, but the present single case cannot prescribe a universal ensemble design.

The model-only d03 Horsepool profiles at 1200 and 1300 UTC (Figure 8) expose the source-specific supplied states and their first-hour evolution. They do not separate early model adjustment from lateral-boundary evolution and do not validate inversion depth.

4.2. The Difference Is a State Offset, Not Demonstrated Cold-Pool Loss

The heat-deficit result narrows the interpretation. GFS had less deficit than NAM under every feedback and footprint definition, but deficit increased in every simulation during the retained window (Figure 7). This sheds doubt that it is erosion of the cold pool (“heat-deficit loss”) that yields a weaker cold pool after spin-up. A longer integration through at least one complete diurnal cycle is required to test whether the products differ in cold-pool erosion or persistence.

The spatial analysis gives a similarly bounded result. The product response was strongest below 500 m AGL and had the expected footprint-mean sign at every hour, but the sign-coherent share fell below two thirds at 1400–1500 UTC. The response is shallow and mean-coherent on d02, not spatially uniform throughout the window.

4.3. Secondary Experiments Narrow Rather than Expand the Claim

The secondary comparisons delimit the driving-product result rather than expanding it into a configuration ranking:

• Figure 4 shows that the driving-product contrasts are furthest from zero on the common temperature, stability, and transport-wind metrics.
• The boundary-layer/surface-layer treatment produced a coherent but smaller response: temperature bias increased while wind absolute error decreased.
• Nest-feedback effects remained small on both default- and fine-terrain branches, providing a five-hour robustness check rather than evidence that feedback is generally unimportant.
• The slope-radiation response was negligible at the station columns, while the supplementary sidewall diagnostic shows a spatially organised response that the station metric missed.
• The controlled terrain brackets show domain- and diagnostic-dependent responses but establish neither a terrain-source ranking nor a drainage or cold-air-volume mechanism.

4.4. Observation and Design Limits Set the Next Experiment

These outcomes motivate better measurements and future experiment design, but they do not support radiation or terrain-source superiority claims. The strongest result is the two-by-two GFS-versus-NAM matrix under one-way and two-way feedback. Two limitations are especially important for the next experiment:

• Native PBL height is scheme dependent. MYJ and YSU need not diagnose the layer top in the same way, so common-height stability and transport remain the more comparable quantities.
• The floor-to-bench temperature contrast is not vertical validation. Independent profiles are required to evaluate inversion depth, layer stability, and modelled heat deficit.

5. Conclusions

The largest repeated response in the controlled archive is the GFS-versus-NAM driving-product contrast on the common d02 evaluation sample. Relative to NAM, GFS is warmer at the evaluated surface sites, about 3.6 K less stable from 50 to 500 m AGL, 0.8–1.0 m s⁻¹ faster in fixed-layer transport wind, and lower in footprint-mean, vertically integrated heat deficit under both feedback modes. The heat-deficit separation is a state offset, not demonstrated faster loss, and the surface metrics do not support a universal product ranking. Nest-feedback effects are much smaller, the PBL/surface-layer response is smaller and metric dependent, and the slope-radiation and fine-terrain responses do not support configuration-superiority claims. The next experiment requires a longer matched matrix, independent vertical observations, a common-support d02–d03 transfer test, and a design that separates initial from lateral-boundary forcing. The present result does not identify the responsible product component, validate inversion depth, establish event-scale persistence, calculate basin export, or predict ozone. Within those bounds, the driving-product treatment produced the largest repeated separation in the lower-basin meteorological state supplied to downstream chemistry during the five-hour window. The numerical summaries supporting the concise Results narrative are retained in Appendix A; Supplementary Tables S1–S2 and Figure S1 retain the bounded secondary comparisons, terrain summary, and slope-aware-radiation response.