Sensitivity of an Idealized Tropical Cyclone to the Configuration of the Global Forecast System–Eddy Diffusivity Mass Flux Planetary Boundary Layer Scheme

The intensity and structure of simulated tropical cyclones (TCs) are known to be sensitive to the planetary boundary layer (PBL) parameterization in numerical weather prediction models. In this paper, we use an idealized version of the Hurricane Weather Research and Forecast system (HWRF) with constant sea-surface temperature (SST) to examine how the configuration of the PBL scheme used in the operational HWRF affects TC intensity change (including rapid intensification) and structure. The configuration changes explored in this study include disabling non-local vertical mixing, changing the coefficients in the stability functions for momentum and heat, and directly modifying the Prandtl number (Pr), which controls the ratio of momentum to heat and moisture exchange in the PBL. Relative to the control simulation, disabling non-local mixing produced a ~15% larger storm that intensified more gradually, while changing the coefficient values used in the stability functions had little effect. Varying Pr within the PBL had the greatest impact, with the largest Pr (~1.6 versus ~0.8) associated with more rapid intensification (~38 versus 29 m s−1 per day) but a 5–10 m s−1 weaker intensity after the initial period of strengthening. This seemingly paradoxical result is likely due to a decrease in the radius of maximum wind (~15 versus 20 km), but smaller enthalpy fluxes, in simulated storms with larger Pr. These results underscore the importance of measuring the vertical eddy diffusivities of momentum, heat, and moisture under high-wind, open-ocean conditions to reduce uncertainty in Pr in the TC PBL.


Introduction
The parameterization of the planetary boundary layer (PBL) in numerical weather prediction models determines how the momentum and enthalpy fluxes from the surface affect the lower atmosphere [1]. It is known that these fluxes play a key role in modulating tropical cyclone (TC) intensity and structure [2][3][4][5], which makes the PBL parameterization critical for achieving accurate TC forecasts. The operational Hurricane Weather Research and Forecast system (HWRF) [6][7][8][9] is NOAA's flagship regional TC prediction system. The atmospheric model includes a variant of the same PBL parameterization-the Global Forecast System-Eddy Diffusivity Mass Flux (GFS-EDMF) scheme-which has been used in the operational Global Forecast System (GFS) since 2015. Readers should note that a newer PBL scheme that predicts turbulent kinetic energy [10] is planned for use in the

Background
The GFS-EDMF PBL consists of both local and non-local vertical mixing parameterizations. The overall diffusive tendency of a prognostic variable X (e.g., zonal and meridional wind speed, potential temperature, or water vapor mixing ratio) [8] is given by ∂X ∂t The first term in brackets in the above equation represents the contribution from the local vertical mixing scheme to the vertical diffusion of X. The local scheme uses K-theory to relate the time tendencies of the prognostic variables to the eddy diffusivity (K X ) and vertical gradient of X. The local scheme is described in more detail below and in previous studies [29,30]. The form of the non-local term F X,NL varies depending on the stability [31]. In weakly unstable conditions, a countergradient term is used to represent the non-local fluxes of heat and momentum [29][30][31]: where c is a coefficient of proportionality assumed to be equal to 6.5 in the GFS-EDMF PBL [29], (w X ) 0 is the surface flux of X, w s is a velocity scale, and h is the PBL height. Under strongly unstable conditions, the GFS-EDMF PBL uses a mass-flux scheme to determine the non-local fluxes [8,31]: where M is the updraft mass flux, X u is the value of X in the updraft, and X is the value of X in the environment. While local mixing is always active, a dimensionless stability parameter (ζ = z 1 /L), where z 1 is the height of the first model vertical level and L is the Obukhov length, is used to determine whether to activate the non-local mixing components, according to the criteria in Table 1. A few other details about the local vertical mixing scheme are relevant to this study. The K m parameterization has received a lot of attention in HWRF and is worth discussing in more depth. K m is calculated from a cubic function that depends on the height above the z α 1 − z h 2 (6) K h and the vertical eddy diffusivity of the water vapor mixing ratio (K q ) are assumed to be equal and are obtained by multiplying K m by the inverse of Pr: K h and K q are ultimately used to calculate the enthalpy fluxes [12]: where SHF is the sensible heat flux, ρ is the air density, c p is the specific heat of air at constant pressure, dθ/dz is the vertical gradient of potential temperature, LHF is the latent heat flux, L v is the latent heat of vaporization, and dq/dz is the vertical gradient of the water vapor mixing ratio. Because Pr is used to determine K h and K q , which in turn are needed to calculate the enthalpy fluxes that sustain the TC, Pr is an important parameter in the GFS-EDMF. When the non-local countergradient (but not the mass-flux) mixing scheme is active, Pr is calculated according to the following equation: where z s is the assumed fractional depth of the surface layer relative to the entire depth of the PBL. Since z s = 0.1 and c = 6.5 in the version of the GFS-EDMF PBL used in HWRF, the entire second term is a constant equal to 0.26. When there is either only local mixing (ζ > 0.2) or the non-local mass-flux scheme is activated (ζ < −0.5), Pr is simply the ratio of the stability function for heat to the stability function for momentum. The stability functions are taken from Dyer and Hicks [32]: where R i is the bulk Richardson number. The first set of stability functions is applied in the unstable regime (R i ≤ 0), while the second set is applied in the stable regime (R i > 0). While Dyer and Hicks [32] use the classical definition of ζ = z/L, the GFS-EDMF PBL uses a different quantity (z s h/L) in place of ζ to calculate the stability functions. The implications of this choice are discussed in Section 4.2.

Methods
The idealized configuration of HWRFv3.9a [8] was used in this study. This version of the model is consistent with the 2017 operational HWRF, but the idealized configuration does not include ocean or wave coupling, vortex initialization, or data assimilation. Instead, the model is initialized with a quiescent sounding that is based on the mean vertical profiles of temperature and humidity from June to September in the Caribbean as reported by Jordan [33]. An idealized vortex, with an initial horizontal wind speed (i.e., intensity) of 20 m s −1 and a radius of maximum wind (RMW) of 90 km, is then superimposed on the base state. The idealized vortex is described in more detail in [34]. There is no land in the model domain, and the sea-surface temperature (SST) is held constant at 302 K (28.85 • C). In all other respects, the idealized HWRFv3.9a configuration and the 2017 operational HWRF are equivalent, and utilize the same 18-6-2 km triply-nested configuration and the same atmospheric physics, which are listed in Table 2 and discussed in more detail in [8]. While the operational HWRF has been upgraded multiple times since 2017, there have not been any changes to the PBL scheme, and most of the remaining changes are related to improving the initialization system (i.e., vortex initialization and data assimilation), which is not utilized here. The simplified initial conditions and the lack of ocean coupling make it easier to isolate the impact of PBL scheme changes on the TC, due to the reduced number of nonlinear interactions relative to the operational system. The first simulation that was conducted utilizes the default configuration of the HWRF PBL that is contained in HWRFv3.9a and is described in the previous section. This control simulation will be referred to as CTRL hereafter. Three modifications to the PBL scheme are then tested in individual experiments. The first change (LOCAL-MIXING) studies the impact of disabling the non-local vertical mixing (i.e., the countergradient and mass-flux components of the vertical mixing) on the TC intensity and structure. The purpose of this experiment is to understand the relative roles of the local and non-local vertical mixing in the TC PBL. We note that some sources of non-local vertical mixing, such as mixing caused by partially-resolved individual updrafts and downdrafts associated with convection, will still be present in the LOCAL-MIXING simulation. Therefore, this experiment only investigates the TC sensitivity to the subgrid-scale (SGS) non-local vertical mixing in the PBL. Relative to CTRL, disabling this SGS non-local vertical mixing in the PBL leads to a smaller Pr, since the second term in Equation (9) is eliminated when countergradient fluxes are not active. Equation (6) suggests that smaller Pr increases K h , which enhances the enthalpy flux from the ocean to the atmosphere for a given near-surface wind speed (Equations (7)-(8)). Therefore, it is reasonable to expect that changes to the HWRF PBL that impact K h might affect the TC intensity and structure.
The second change (BUSINGER-PHI) replaces the stability functions for momentum and heat from Dyer and Hicks [32] (Equations (10)-(13)), with an alternative form proposed by Businger et al. [35] (Equations (14)-(17)): ϕ m and ϕ h are used to calculate K m , K h , and Pr via Equations (5), (6) and (9), respectively. The effect of this change is to decrease Pr slightly, but to a much lesser extent than in LOCAL-MIXING, which we will show later. In a final experiment (Pr = 1.58), the Pr within the PBL is set to a constant value of 1.58, which is much larger than that suggested by either the Dyer or Businger relations, to bookend the results from the other simulations. This specific Pr value was chosen because observations of K m and K h in the TC PBL collected during the Coupled Boundary Layer Air-Sea Transfer (CBLAST) experiment imply a mean Pr of 1.58 [12]. Above the PBL, the formulation for Pr was not changed. Because there is substantial scatter in the individual values of Pr = K m /K h from CBLAST, this experiment is undertaken as a sensitivity test, and is not meant to imply that a Pr of 1.58 is a reasonable choice throughout the TC PBL. This point is discussed in more detail in Section 4.4.
While the above experiments may appear to be unrelated, we will show later that each of their results can be explained by understanding their impact on Pr and K h . Table 3 summarizes the control and experiment simulations. In each experiment, the GFDL vortex tracker [36] was used in the inner nest (vor-tex_tracker option 7 in the namelist) to track the storm center position, maximum 10-m wind speed, and minimum sea-level pressure. The parameters used to estimate the storm center location were the (1) relative vorticity at 10 m, 850 hPa, and 700 hPa; (2) horizontal wind speed at 10 m, 850 hPa, and 700 hPa; (3) mean sea-level pressure; and (4) geopotential height at 850 hPa and 700 hPa. The version of the tracker used here ran concurrently with the HWRF forecast, providing updated storm information every 6 min. The inline tracker is identical to the external tracker that runs as a postprocessing step in the operational HWRF, with the exception that the internal tracker does not attempt to detect TC dissipation. Since dissipation does not occur in any of our simulations, we expect our results to be comparable to those produced by the external tracker.
While the GFDL vortex tracker is capable of providing wind radii information, we chose to calculate our own radii of 34-kt (tropical-storm force), 64-kt (hurricane-force), and maximum winds (1 kt = 0.514 m s −1 ). For this purpose, we transformed the threedimensional model output from the inner nest to cylindrical coordinates (relative to the TC center) and calculated the wind radii as a function of height from the azimuthally-averaged wind field. The GFDL vortex tracker, in contrast, only provides the wind radii at 10 m. In addition, the tracker calculates the maximum extent of each wind radius, rather than the values from the azimuthally-averaged wind field that we computed. Considering the wind radii of the azimuthally-averaged wind field reduces noise in the time series, making it easier to visualize differences in the storm size between the various experiments.

Characteristics of the CTRL Simulation
Shortly after initialization, the TC in CTRL has an intensity of~15 m s −1 ( Figure 1a) and a RMW of 90 km (Figure 1c). During the first day of the simulation, a period of rapid intensification [37] ensues, with the RMW contracting to~20 km and the storm intensity increasing to~53 m s −1 by t = 24 h. A brief period of weakening occurs from t = 40-60 h before the storm reaches a second, stronger peak intensity of~70 m s −1 around t = 100 h. Some fluctuations in intensity occur during the final day of the simulation (t = 100-120 h). There is generally close correspondence between the trends in maximum wind speed and minimum central pressure ( Figure 1b). The storm size increases as the TC matures, with the radius of 34-kt winds (dotted black line in Figure 1c) approximately doubling from t = 20 to 120 h. ter) and calculated the wind radii as a function of height from the azimuthally-averaged wind field. The GFDL vortex tracker, in contrast, only provides the wind radii at 10 m. In addition, the tracker calculates the maximum extent of each wind radius, rather than the values from the azimuthally-averaged wind field that we computed. Considering the wind radii of the azimuthally-averaged wind field reduces noise in the time series, making it easier to visualize differences in the storm size between the various experiments.

Characteristics of the CTRL Simulation
Shortly after initialization, the TC in CTRL has an intensity of ~15 m s −1 ( Figure 1a) and a RMW of 90 km (Figure 1c). During the first day of the simulation, a period of rapid intensification [37] ensues, with the RMW contracting to ~20 km and the storm intensity increasing to ~53 m s −1 by t = 24 h. A brief period of weakening occurs from t = 40-60 h before the storm reaches a second, stronger peak intensity of ~70 m s −1 around t = 100 h. Some fluctuations in intensity occur during the final day of the simulation (t = 100-120 h). There is generally close correspondence between the trends in maximum wind speed and minimum central pressure ( Figure 1b). The storm size increases as the TC matures, with the radius of 34-kt winds (dotted black line in Figure 1c) approximately doubling from t = 20 to 120 h. Because the value of ζ = z 1 /L determines whether and what kind of non-local mixing is active within a particular vertical column of the model simulation (Table 1), it is useful to examine the spatial distribution of z 1 /L within the TC to understand the behavior of the GFS-EDMF scheme there. Figure 2 shows the value of z 1 /L on the innermost HWRF domain for three different forecast hours in CTRL during the mature stage of the TC. The 34-, 50-, and 64-kt wind contours (black lines in Figure 2) are also shown to demonstrate how the model values of z 1 /L correspond to the TC structure. In general, stability increases with decreasing distance from the TC center, which is consistent with the assumption that the TC PBL is near neutral (z 1 /L~0). The value of z 1 /L is less than 0.2 everywhere in the domain, indicating that countergradient fluxes are active both inside the TC and within  (Table 1). However, there are only small, isolated areas where z 1 /L is less than −0.5 and the mass flux scheme is active (white areas in Figure 2). The regions that have z 1 /L near or less than −0.5 are all associated with rainbands. This result suggests that the mass flux scheme probably does not play a substantial role in modulating TC structure or intensity.
to examine the spatial distribution of z1/L within the TC to understand the behavior of the GFS-EDMF scheme there. Figure 2 shows the value of z1/L on the innermost HWRF domain for three different forecast hours in CTRL during the mature stage of the TC. The 34-, 50-, and 64-kt wind contours (black lines in Figure 2) are also shown to demonstrate how the model values of z1/L correspond to the TC structure. In general, stability increases with decreasing distance from the TC center, which is consistent with the assumption that the TC PBL is near neutral (z1/L~0). The value of z1/L is less than 0.2 everywhere in the domain, indicating that countergradient fluxes are active both inside the TC and within the environment (Table 1). However, there are only small, isolated areas where z1/L is less than −0.5 and the mass flux scheme is active (white areas in Figure 2). The regions that have z1/L near or less than −0.5 are all associated with rainbands. This result suggests that the mass flux scheme probably does not play a substantial role in modulating TC structure or intensity. The maintenance of positive temperature and moisture contrasts between the ocean and the atmosphere (sea minus air) is essential for sustaining the enthalpy fluxes that drive the TC heat engine. Here, we examine the difference in the air temperature (dT; Figure 3a) and water vapor mixing ratios (dq; Figure 4a) between the ocean surface and the first model vertical level (approximately 10-m height), since those are the quantities and vertical levels that the GFDL surface layer scheme uses to calculate the fluxes from the ocean to the atmosphere. To calculate the water vapor mixing ratio at the ocean surface, we assume that the SST is equal to the air temperature and that the air is saturated (i.e., 100% relative humidity). Within the 50-kt wind radius, dT and dq range from 3 to 4 K ( Figure 3a) and 6-7.5 g kg −1 (Figure 4a), respectively, at t = 120 h in CTRL. The model values are about 50% larger than buoy observations of dT and dq at 10-m height in hurricanes [38], which is consistent with the overly-unstable HWRF PBL noted in previous studies [23,39]. The shape of the observed and simulated radial distributions of dT and dq, however, are qualitatively similar. These values of dT and dq contribute to total (sensible + latent) enthalpy fluxes of up to 1850 W m −2 within the 50-kt wind radius (Figure 5a), which we refer to as the TC inner core in the remainder of this study. The maintenance of positive temperature and moisture contrasts between the ocean and the atmosphere (sea minus air) is essential for sustaining the enthalpy fluxes that drive the TC heat engine. Here, we examine the difference in the air temperature (dT; Figure 3a) and water vapor mixing ratios (dq; Figure 4a) between the ocean surface and the first model vertical level (approximately 10-m height), since those are the quantities and vertical levels that the GFDL surface layer scheme uses to calculate the fluxes from the ocean to the atmosphere. To calculate the water vapor mixing ratio at the ocean surface, we assume that the SST is equal to the air temperature and that the air is saturated (i.e., 100% relative humidity). Within the 50-kt wind radius, dT and dq range from 3 to 4 K ( Figure 3a) and 6-7.5 g kg −1 (Figure 4a), respectively, at t = 120 h in CTRL. The model values are about 50% larger than buoy observations of dT and dq at 10-m height in hurricanes [38], which is consistent with the overly-unstable HWRF PBL noted in previous studies [23,39]. The shape of the observed and simulated radial distributions of dT and dq, however, are qualitatively similar. These values of dT and dq contribute to total (sensible + latent) enthalpy fluxes of up to 1850 W m −2 within the 50-kt wind radius (Figure 5a), which we refer to as the TC inner core in the remainder of this study.

Non-local Mixing and Simulated Tropical Cyclones
To determine whether TCs are sensitive to non-local mixing in the GFS-EDMF PBL, we performed an experiment (LOCAL-MIXING) in which the mass-flux scheme was deactivated and the countergradient term was set to zero throughout the model domain. This leaves only the local mixing scheme active. Figure 1a demonstrates that when nonlocal mixing is disabled, the TC intensifies more gradually, reaching a weaker initial peak intensity (58 m s −1 ) than in CTRL (63 m s −1 ). However, the storm in LOCAL-MIXING eventually obtains the same peak wind speed as the CTRL storm, and maintains this wind speed, along with a smaller minimum central pressure (Figure 1b), for a longer portion of the t = 100-120 h period. An additional difference between the CTRL and LOCAL-MIXING storms is that the various wind radii are larger in LOCAL-MIXING, with a~15% larger 34-kt wind radius (Figure 1c).  .   .

Non-local Mixing and Simulated Tropical Cyclones
To determine whether TCs are sensitive to non-local mixing in the GFS-EDMF PBL, we performed an experiment (LOCAL-MIXING) in which the mass-flux scheme was deactivated and the countergradient term was set to zero throughout the model domain. This leaves only the local mixing scheme active. Figure 1a demonstrates that when nonlocal mixing is disabled, the TC intensifies more gradually, reaching a weaker initial peak intensity (58 m s −1 ) than in CTRL (63 m s −1 ). However, the storm in LOCAL-MIXING eventually obtains the same peak wind speed as the CTRL storm, and maintains this wind speed, along with a smaller minimum central pressure (Figure 1b), for a longer portion of the t = 100-120 h period. An additional difference between the CTRL and LOCAL-MIX-ING storms is that the various wind radii are larger in LOCAL-MIXING, with a ~15% larger 34-kt wind radius (Figure 1c).
Relative to the CTRL storm, the total enthalpy flux in the LOCAL-MIXING storm is 15-50% larger within the 50-kt wind radius at t = 120 h (Figure 5b). At t = 120 h, the CTRL and LOCAL-MIXING storms have the same intensity (Figure 1a), which indicates that the difference in the enthalpy fluxes is mainly driven by differences in dT and dq rather than wind speed. Indeed, dT ( Figure 3b) and dq (Figure 4b) are roughly 15-40% larger within the 50-kt wind radius in the LOCAL-MIXING simulation. The larger sea-air contrasts and greater enthalpy fluxes in LOCAL-MIXING are explained by the reduction in Pr that occurs in this experiment (Figure 6b) relative to CTRL (Figure 6a). The smaller Pr in LOCAL-MIXING is due to disabling the countergradient fluxes, which eliminates the second term in the Pr calculation (Equation (9)) and reduces Pr by 0.26 when −0.5 ≤ ≤ 0.2. The smaller Pr leads to increased Kh (Equation (6)), Kq (Equation (6)), and enthalpy fluxes in the PBL (Equations (7)-(8)). Recall that the dimensionless ratio used to calculate the stability functions in the GFS-EDMF is ℎ/ , not . The effect of this choice is to reduce Pr below the value implied by the Dyer relation (black line in Figure 6b) for a given value of . In CTRL, the addition of the second term in Equation (9) boosts Pr close to the values suggested by Relative to the CTRL storm, the total enthalpy flux in the LOCAL-MIXING storm is 15-50% larger within the 50-kt wind radius at t = 120 h (Figure 5b). At t = 120 h, the CTRL and LOCAL-MIXING storms have the same intensity (Figure 1a), which indicates that the difference in the enthalpy fluxes is mainly driven by differences in dT and dq rather than wind speed. Indeed, dT ( Figure 3b) and dq (Figure 4b) are roughly 15-40% larger within the 50-kt wind radius in the LOCAL-MIXING simulation. The larger sea-air contrasts and greater enthalpy fluxes in LOCAL-MIXING are explained by the reduction in Pr that occurs in this experiment (Figure 6b) relative to CTRL (Figure 6a). The smaller Pr in LOCAL-MIXING is due to disabling the countergradient fluxes, which eliminates the second term in the Pr calculation (Equation (9)) and reduces Pr by 0.26 when −0.5 ≤ ζ ≤ 0.2. The smaller Pr leads to increased K h (Equation (6)), K q (Equation (6)), and enthalpy fluxes in the PBL (Equations (7)-(8)). Recall that the dimensionless ratio used to calculate the stability functions in the GFS-EDMF is z s h/L, not ζ. The effect of this choice is to reduce Pr below the value implied by the Dyer relation (black line in Figure 6b) for a given value of ζ. In CTRL, the addition of the second term in Equation (9) boosts Pr close to the values suggested by the Dyer relation, but only when countergradient fluxes are active (Figure 6a). The relationship between Pr, K h , and the enthalpy fluxes will be discussed in greater detail in Section 4.4.
Atmosphere 2021, 12, x FOR PEER REVIEW 11 of 18 the Dyer relation, but only when countergradient fluxes are active (Figure 6a). The relationship between Pr, Kh, and the enthalpy fluxes will be discussed in greater detail in Section 4.4.

Stability Functions and Simulated Tropical Cyclones
Within the GFS-EDMF PBL, the stability functions for heat and momentum are taken from Dyer (Section 2; Equations. (10)-(13)). These stability functions are used within the PBL to calculate Pr from Equation (9). In a third experiment (BUSINGER-PHI), Dyer's coefficients were replaced with Businger's (Equations (14)- (17)) in the GFS-EDMF PBL to determine whether this particular choice of stability functions impacts TC intensity and structure. BUSINGER-PHI produces a storm that has a slightly weaker peak intensity, with a maximum wind that is ~5 m s −1 weaker than CTRL for most forecast lead times from 80 to 120 hours (Figure 1a). The minimum central pressure is also about 5 hPa larger during most of the simulation after the first 36 h (Figure 1b). In general, the storm created in BUSINGER-PHI also does not intensify as rapidly as the storm in CTRL and takes almost twice as long to exceed an intensity of 60 m s −1 . The CTRL and BUSINGER-PHI storms have fairly similar wind radii (Figure 1c). Overall, dT (Figure 3c), dq (Figure 4c), and the total enthalpy flux (Figure 5c) do not show a consistent change relative to CTRL in the TC inner core. Unlike in the LOCAL-MIXING experiment, Pr is only slightly smaller in the BUSINGER-PHI experiment than in CTRL (by about 0.1; Figure 6c), which contributes to the similar sea-air contrasts and enthalpy fluxes in the latter two simulations. Figure 6c also indicates that the mass-flux scheme is more active in the TC environment in

Stability Functions and Simulated Tropical Cyclones
Within the GFS-EDMF PBL, the stability functions for heat and momentum are taken from Dyer (Section 2; Equations. (10)-(13)). These stability functions are used within the PBL to calculate Pr from Equation (9). In a third experiment (BUSINGER-PHI), Dyer's coefficients were replaced with Businger's (Equations (14)- (17)) in the GFS-EDMF PBL to determine whether this particular choice of stability functions impacts TC intensity and structure. BUSINGER-PHI produces a storm that has a slightly weaker peak intensity, with a maximum wind that is~5 m s −1 weaker than CTRL for most forecast lead times from 80 to 120 hours (Figure 1a). The minimum central pressure is also about 5 hPa larger during most of the simulation after the first 36 h (Figure 1b). In general, the storm created in BUSINGER-PHI also does not intensify as rapidly as the storm in CTRL and takes almost twice as long to exceed an intensity of 60 m s −1 . The CTRL and BUSINGER-PHI storms have fairly similar wind radii (Figure 1c). Overall, dT (Figure 3c), dq (Figure 4c), and the total enthalpy flux (Figure 5c) do not show a consistent change relative to CTRL in the TC inner core. Unlike in the LOCAL-MIXING experiment, Pr is only slightly smaller in the BUSINGER-PHI experiment than in CTRL (by about 0.1; Figure 6c), which contributes to the similar sea-air contrasts and enthalpy fluxes in the latter two simulations. Figure 6c also indicates that the mass-flux scheme is more active in the TC environment in BUSINGER-PHI than in CTRL, since there are more points with weaker than tropical-storm-force winds that have ζ < −0.5.

Tying the Results Together: Prandtl Number and Tropical Cyclones
The CBLAST experiment derived Pr from observations of K m and K h within the TC PBL in wind speeds that ranged from 18 to 30 m s −1 [12]. Their observations suggest that Pr values below one, such as those in the LOCAL-MIXING experiment (Figure 6b), are too small for the hurricane environment. In fact, their K m and K h observations imply a mean Pr of 1.58 (Figure 9 in Zhang and Drennan [12]), which is much larger than the values in CTRL or in the experiments conducted so far. However, the applicability of these data to our model simulations is questionable, since the observations were made in multiple TCs at heights from 60 to 800 m in conditions that varied from weakly unstable to stable, as determined from the vertical gradient of potential temperature measured during the flights. In addition, there is substantial scatter in each individual value of Pr = K m /K h (their Figure 9), which suggests that the mean Pr value may not be particularly informative. For these reasons, we do not wish to imply that Pr = 1.58 is a realistic value to use throughout the TC PBL. Rather, our purpose is to determine the impact of using a larger Pr value on the structure and evolution of the TC. To consider this effect, a fourth experiment (Pr = 1.58) was conducted in which Pr was set to a fixed value of 1.58 within the PBL. Above the PBL, Pr was not changed.
The simulated maximum wind speed, minimum central pressure, and RMW from this experiment are shown by the green curves in Figure 1. Compared to the other three experiments, the TC in Pr = 1.58 undergoes the most rapid rate of deepening, reaching a sub-950-hPa pressure approximately ten hours earlier than the storms in the other experiments (Figure 1b). This rapid deepening is also reflected in the maximum wind. Notably, however, the minimum central pressure and maximum winds level off at values corresponding to a weaker storm than in the other experiments. In addition, the storm size, especially the 34-kt wind radius, is much smaller than in the control (100 km versus 125 km at t = 120 h).
How are the smaller storm size, greater propensity for rapid intensification, and weaker peak intensity in the Pr = 1.58 storm related? The comparison of the enthalpy fluxes in the four simulations suggests an explanation. The inner-core enthalpy fluxes at t = 120 h (Figure 5d) are up to 45% weaker in the Pr = 1.58 simulation than in CTRL. Even in the environment, the fluxes are notably weaker. While enthalpy fluxes are generally strongly correlated with wind speed, differences in wind speed do not explain why the enthalpy fluxes differ between these particular simulations, since at t = 120 h the storms have similar maximum wind speeds (Figure 1a). Instead, the air-sea temperature ( Figure 3) and moisture contrasts ( Figure 4) are substantially reduced in the Pr = 1.58 storm, which in turn reduces the enthalpy fluxes within both the environment and the TC inner core. In addition, because K m is unchanged, the larger Pr (=K m /K h ) implies smaller K h for the same wind speed (Figure 7), which further reduces the enthalpy fluxes in the PBL via Equations (7)- (8). Regarding the smaller storm size in Pr = 1.58, controls on TC structure are not well understood, but a few studies have posited a relationship between storm size and the magnitude and aerial extent of the surface enthalpy fluxes [40,41]. If this relationship is valid, it may explain why the Pr = 1.58 storm is smaller than the storm in CTRL.
Regardless of the cause for the more compact TC in Pr = 1.58, smaller storms tend to intensify more rapidly [42][43][44][45], and the storm in Pr = 1.58 strengthens faster than the TCs in any of the other simulations (Figure 1). When the initial period of strengthening begins in Pr = 1.58 at around t = 10 h, the storm is weak with a peak wind of~15 m s −1 , and it is presumed that large enthalpy fluxes are not needed to intensify such a weak TC. The RMW in Pr = 1.58, however, is already beginning to decrease below the RMW of the CTRL storm, and this trend continues throughout the RI period. Near the end of the RI period at t = 24 h, the maximum winds of the Pr = 1.58 storm have already exceeded 50 m s −1 . As the wind speed approaches the maximum intensity, the reduced enthalpy fluxes in the Pr = 1.58 storm play a larger role, and it reaches an earlier but weaker initial peak intensity (~57 m s −1 ) than the storm in CTRL (63 m s −1 ).
Atmosphere 2021, 12, x FOR PEER REVIEW 13 of 18 = 1.58 storm play a larger role, and it reaches an earlier but weaker initial peak intensity (~57 m s −1 ) than the storm in CTRL (63 m s −1 ). The Pr = 1.58 simulation is not the only experiment in which the relationship between Pr, enthalpy fluxes, storm size, intensification rate, and peak intensity is evident. Instead, the results from the LOCAL-MIXING experiment can also be explained through this mechanism. Figure 8 demonstrates that near the RMW within the PBL, the storm in the LOCAL-MIXING experiment has a Pr value of ~0.8, respectively, compared to 1.1 in CTRL (Figure 8). Within the TC inner core, LOCAL-MIXING also has larger sea-air contrasts (Figures 3 and 4), Kh values at a given wind speed (Figure 7), and enthalpy fluxes ( Figure  5) than CTRL, possibly leading to a larger storm size (Figure 1c). The LOCAL-MIXING TC intensifies more slowly than the CTRL and Pr = 1.58 storms during the first 36 h due to its larger size, but ultimately reaches a peak intensity that is comparable to or greater than the CTRL and Pr = 1.58 storms due to its larger enthalpy fluxes. In contrast, the changes to the stability functions in the BUSINGER-PHI experiment have a relatively minor impact on Pr (Figures 6c and 8c), which explains why the CTRL and BUSINGER-PHI TCs are similar. The Pr = 1.58 simulation is not the only experiment in which the relationship between Pr, enthalpy fluxes, storm size, intensification rate, and peak intensity is evident. Instead, the results from the LOCAL-MIXING experiment can also be explained through this mechanism. Figure 8 demonstrates that near the RMW within the PBL, the storm in the LOCAL-MIXING experiment has a Pr value of~0.8, respectively, compared to 1.1 in CTRL (Figure 8). Within the TC inner core, LOCAL-MIXING also has larger sea-air contrasts (Figures 3 and 4), K h values at a given wind speed (Figure 7), and enthalpy fluxes ( Figure 5) than CTRL, possibly leading to a larger storm size (Figure 1c). The LOCAL-MIXING TC intensifies more slowly than the CTRL and Pr = 1.58 storms during the first 36 h due to its larger size, but ultimately reaches a peak intensity that is comparable to or greater than the CTRL and Pr = 1.58 storms due to its larger enthalpy fluxes. In contrast, the changes to the stability functions in the BUSINGER-PHI experiment have a relatively minor impact on Pr (Figures 6c and 8c), which explains why the CTRL and BUSINGER-PHI TCs are similar. Atmosphere 2021, 12, x FOR PEER REVIEW 14 of 18 Although the impact of Pr on TC intensity change is evident, the most striking difference is the impact of Pr on the size of the tropical storm-force wind radius at 10-m height. At t = 120 h, the 34-kt wind radius is 100 km in the experiment with the largest Pr (Pr = 1.58), but 150 km in the experiment with the smallest Pr (0.8, LOCAL-MIXING), a remarkable 40% difference. Recent (2018-2020) versions of the operational HWRF exhibit a negative 34-kt wind radius bias, a near-neutral 50-kt wind radius bias, and positive 64kt wind radius and RMW biases [46], so it is not clear that adjusting Pr in the operational system should be considered, since all of the wind radii will be affected. However, newer TC forecasting systems under development at the National Centers for Environmental Prediction that also use the GFS-EDMF PBL, such as the Hurricane Analysis and Forecast System, currently exhibit positive wind radii biases [47]. Adjustments to Pr in these models could be considered if they are done in a physically-consistent manner and do not degrade the TC track and intensity forecasts.

Summary and Conclusions
In this paper, we used an idealized version of HWRF to study the sensitivity of TC intensity and structure to several settings in the GFS-EDMF PBL, which is used in the operational HWRF. The model was initialized with an idealized vortex that was superimposed on a quiescent base state and ran for five days. The SST was fixed at 302 K (28.85 °C) Although the impact of Pr on TC intensity change is evident, the most striking difference is the impact of Pr on the size of the tropical storm-force wind radius at 10-m height. At t = 120 h, the 34-kt wind radius is 100 km in the experiment with the largest Pr (Pr = 1.58), but 150 km in the experiment with the smallest Pr (0.8, LOCAL-MIXING), a remarkable 40% difference. Recent (2018-2020) versions of the operational HWRF exhibit a negative 34-kt wind radius bias, a near-neutral 50-kt wind radius bias, and positive 64-kt wind radius and RMW biases [46], so it is not clear that adjusting Pr in the operational system should be considered, since all of the wind radii will be affected. However, newer TC forecasting systems under development at the National Centers for Environmental Prediction that also use the GFS-EDMF PBL, such as the Hurricane Analysis and Forecast System, currently exhibit positive wind radii biases [47]. Adjustments to Pr in these models could be considered if they are done in a physically-consistent manner and do not degrade the TC track and intensity forecasts.

Summary and Conclusions
In this paper, we used an idealized version of HWRF to study the sensitivity of TC intensity and structure to several settings in the GFS-EDMF PBL, which is used in the operational HWRF. The model was initialized with an idealized vortex that was superimposed on a quiescent base state and ran for five days. The SST was fixed at 302 K (28.85 • C) throughout the simulation. The experiments included (1) disabling the non-local vertical mixing parameterizations in the PBL (LOCAL-MIXING), (2) replacing the coefficient values in the stability functions for momentum and heat with those proposed by Businger [35] (BUSINGER-PHI), and (3) using a larger value of the Prandtl number within the PBL (Pr = 1.58), consistent with TC PBL observations from the CBLAST experiment. The results from these experiments were compared to each other and to a control simulation (CTRL) with no changes.
The CTRL simulation produced a TC that underwent a period of rapid intensification (27 m s −1 in 24 h) during the first day, peaking at 70 m s −1 at t = 100 h before weakening slowly thereafter. The storm in the LOCAL-MIXING experiment intensified more slowly than the CTRL storm, but ultimately reached a similar peak intensity and maintained this intensity for a longer period of time. The modifications made in the LOCAL-MIXING experiment (i.e., disabling the non-local vertical mixing parameterizations in the GFS-EDMF) decreased Pr, which resulted in sea-air temperature and moisture contrasts that were 15-40% larger than in the CTRL storm and contributed to a 15-50% larger total enthalpy flux under similar wind conditions. The wind radii in the LOCAL-MIXING storm were~15% larger than those in the CTRL simulation, possibly due to the larger enthalpy fluxes. Because the storm was larger, it intensified more slowly than the TC in CTRL.
In contrast, the larger Pr in the Pr = 1.58 storm reduced the temperature and moisture contrast between the ocean surface and the lowest model level (~10 m) by 30-40%, which decreased the enthalpy fluxes and resulted in a slightly weaker peak intensity (~65 m s −1 ). All of the wind radii in the Pr = 1.58 storm were substantially smaller than those of the CTRL or LOCAL-MIXING storms, especially the 34-kt wind radius, which was 25% and 33% smaller than in the CTRL and LOCAL-MIXING storms, respectively. Despite the weaker peak intensity, the smaller size of the Pr = 1.58 storm allowed it to intensify even more rapidly in the first 24 h than the CTRL or LOCAL-MIXING storms. Unlike the LOCAL-MIXING or Pr = 1.58 experiments, modifying the stability functions in the BUSINGER-PHI experiment had a relatively small impact on Pr and the TC. The most notable difference in the BUSINGER-PHI experiment was that the storm was about 5 m s −1 weaker than the CTRL storm, despite similar sea-air contrasts.
The above results suggest that the impact to Pr should be considered when future changes are made to the GFS-EDMF PBL. In the GFS-EDMF, when Pr is increased, the resulting TC: • Has a smaller K h and K q for a given wind speed; • Exhibits reduced enthalpy fluxes within its inner-core, even at similar wind speeds; • Has smaller radii of maximum, 64-, and 34-kt winds, possibly due to the weaker enthalpy fluxes, with the largest impact on the 34-kt wind radius; • Is more likely to undergo an initial period of rapid intensification, due to its smaller size; • Reaches a weaker peak intensity due to the smaller enthalpy fluxes in the inner core.
Due to the way that the eddy diffusivities of momentum, heat, and moisture are calculated in the GFS-EDMF PBL (Equations (5)-(6)), the eddy diffusivity of momentum for a particular wind speed is unaffected by modifications to Pr, since Pr is only used to obtain the eddy diffusivities of heat and moisture from K m . While K m plays an important role in determining TC intensity and the potential for rapid intensification, it is far from the only important variable in the PBL scheme, as this study demonstrates.
Future work should focus on collecting additional observations of K m and K h in the TC PBL that can be used to determine Pr, possibly using Uncrewed Aircraft Systems equipped with turbulence probes. These additional measurements are needed to further explore the relationship between stability and Pr in the TC PBL, which may differ from that in the more readily-studied PBL over the land, as suggested by the CBLAST results. If subsequent measurements confirm that Pr in the TC PBL is unique, a Pr parameterization that accounts for this behavior could be developed. In addition to stability, such a parameterization might consider whether the vertical column is over land or over water and/or include a wind speed dependence to estimate Pr accurately across a variety of phenomena, including TCs. In addition, examining the TC response across a wider range of Pr, particularly at smaller values than those examined here (i.e., Pr < 0.8), would help to quantify the TC sensitivity to Pr more completely.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to their large size.