Introduction

One of the greatest challenges in numerical weather prediction (NWP) is how to appropriately treat the sub-grid-scale (SGS) transport induced by the unresolved turbulent eddies. The problem is intrinsic to numerical models that utilize discretized grids to simulate the continuous atmospheric flow. Because of the nonlinear interaction between the resolved and unresolved motions, inclusion of SGS processes in the grid-box-mean governing equations results in new second-order terms, causing the otherwise closed forecasting system described by the physical laws to be no longer closed. This is the notorious closure problem for any dynamic models for simulating turbulent flow. To close the system, the second-order moments must be parameterized based on the model-resolved variables. The nature of the problem may be understood by Kolmogorov’s spectrum of turbulent kinetic energy (TKE)1. The turbulent energy peaks at large eddies with scales of a few hundred meters2. These energy-containing eddies are generated by the instabilities of the mean flow and obtain energy directly from the mean flow, and thus, they are fundamentally anisotropic. With the decrease of scales, eddies become less flow dependent and contain less energy. On the energy spectrum, there exists an intermediate range of scales, known as the Kolmogorov inertial subrange1, where the net incoming energy from larger eddies is in equilibrium with the net energy cascading down to smaller eddies. Eddies with sizes smaller than the inertial subrange are commonly considered to be isotropic.

Because of the distinct characteristics between anisotropic large and isotropic small eddies, parameterization of SGS turbulence is a scale-sensitive problem. At one extreme where the model grid spacing falls in or below the inertial subrange, the energy-containing eddies are explicitly resolved. The remaining unresolved isotropic small eddies need to be treated three-dimensionally (3D) to account for their induced directionally independent mixing. This type of high-resolution simulations is known as the large eddy simulations (LESs). Most common 3D turbulence schemes used in LESs include the 3D Smagorinsky–Lilly scheme3 and the 3D Deardorff TKE scheme4. A short review of classic 3D schemes for LESs is provided in Section 4: Methods.

At the other extreme as the case of large-scale and mesoscale (LSMS) models where model resolution is too coarse to resolve energy-containing eddies, parameterization needs to account for the directionally different turbulent mixing induced by anisotropic large eddies. The current method is to treat horizontal and vertical turbulent mixing separately in the model’s dynamic solver and a standalone one-dimensional (1D) module outside the dynamic core, often known as the planetary boundary layer (PBL) scheme. There are scientific and technical reasons for such separate treatments of horizontal and vertical turbulent mixing. Technically, since most model physics, such as microphysics and radiation, are parameterized in a 1D vertical column framework, it is natural to parameterize vertical turbulent mixing in the same 1D framework so that different model physics can be bundled together in a so-called model physics package outside the model dynamic core. However, since this 1D column framework only allows for vertical gradient calculation, horizonal turbulent mixing must be treated within the dynamic solver where the required horizontal gradients can be calculated. This technical separation of vertical and horizontal turbulent mixing parameterization under a model’s infrastructure is supported by scientific considerations. First, since the vertical gradients of grid-box mean variables are normally much larger than the corresponding horizontal gradients at low resolutions, it is scientifically sound to separate the vertical and horizontal turbulent mixing parameterization. Second, there are unique vertical turbulent processes associated with shallow clouds (stratocumulus and fair-weather cumulus) and non-local mixing induced by large turbulent eddies that require special treatment. These processes are relatively easy to handle in a 1D column framework.

In between the LES and LSMS extremes is the so-called gray zone where energy-containing eddies are partially resolved. Classic 3D turbulence schemes designed for LESs cannot be extended to the gray zone since the unresolved eddies are no longer isotropic. To date, all full-physics dynamic models for NWP retain the separated parameterization strategy for treating vertical and horizontal turbulent mixing. But the methodology becomes problematic since in the gray zone the 3D anisotropic turbulent eddies tend to generate inter-connected horizontal and vertical mixing that cannot be artificially separated.

An integrated parameterization for inter-connected anisotropic 3D turbulent mixing is particularly important to numerical simulations of phenomena that consist of both deep and shallow convections. One example is the tropical cyclones (TCs). Figure 1 shows the distinct vertical thermodynamic structure in different regimes of TC inner core. The thermodynamic structure in the TC eye and the moat in between the eyewall and rainbands shares similar characteristics to that of typical fair-weather conditions in that the turbulent boundary layer is separated from the free atmosphere above by a capping inversion layer. In the deep convective eyewall and rainbands, however, no such interface exists to separate the turbulence generated by the boundary layer processes and cloud processes aloft. The boundary layer there becomes ill-defined as the air is pulled up into the active convection5,6,7,8,9. Large TKE extending from the surface to upper troposphere in the eyewall and rainbands has been shown by numerical simulations10,11 and confirmed by the airborne Doppler radar observations12,13. The thermodynamic structural change in the transition from deep to shallow convection regimes implies that the 3D turbulent mixing is important there. Moreover, both observations and LESs show that turbulent eddies experience large lateral gradients across the eyewall, rainbands, and the moat in between (Supplementary Note 1 and Supplementary Fig. 1), resulting in comparable and inter-connected horizontal and vertical turbulent mixing14,15. Such an inter-connected 3D turbulent mixing poses a challenge for the separated vertical and horizontal turbulence parameterizations currently used in operational forecasting models. Finally, as model resolution increases, the 3D shear production of TKE becomes important (Supplementary Note 2 and Supplementary Fig. 2), but it is not considered in the separated turbulent mixing parameterizations.

Fig. 1: Potential temperature (θ) profiles in different regimes of TC inner core.
figure 1

a–d Profiles of potential temperature (\(\theta\)) from 20 individual dropsondes collected within the radii of 150 km from the storm center of Hurricanes Earl (2010), Harvey (2017), Irma (2017), Dorian (2019), and Ida (2021). The soundings are classified into storm eye, eyewall, moat, and rainbands based on the wind speed, relative humidity, and radial distance to the storm center. Note that for Harvey, Irma, Dorian, and Ida the reconnaissance flight level is about 3 km, therefore, soundings above that are not available.

Full size image

There have been efforts to develop 3D turbulence parameterization for the gray zone. These include extending the 3D Smagorinsky–Lilly model beyond the LES limit16, using a simple weighting function to blend a conventional 1D non-local PBL scheme with the 3D Smagorinsky–Lilly model17, merging the LES and LSMS extremes using a scale-aware (SA) mixing length formulation18, and developing a 3D-SA TKE scheme by blending the 3D Deardorff TKE model with a 1D PBL scheme including the non-local mixing effects19. However, many issues regarding 3D turbulence parameterization applicable to the gray zone remain to be addressed. First, previous studies only focused on the dry convective boundary layer (CBL) in horizontally homogenous conditions. It has never been explored how anisotropic 3D turbulent mixing can be realistically parameterized in different weather conditions including weather extremes, such as TCs where different convective regimes possess distinct thermodynamic and turbulence structures. Second, no 3D-SA turbulence scheme has ever been developed for operational NWP models and tested in real-time or near real-time operational forecasts. Finally, one of NOAA’s ongoing strategic efforts is to transition from maintaining multiple, separate weather prediction systems to a single Unified Forecast System (UFS, https://www.ufs.epic.noaa.gov/). UFS with the Finite-Volume Cubed-Sphere Dynamical Core (FV3)20 is a paradigm shift that encourages and promotes a collaborative, community-based development approach to accelerate research and improve the accuracy of NWP across scales. This motivates us to develop a 3D-SA TKE scheme compatible with the FV3 dynamic core and implement it in UFS’s Hurricane Analysis and Forecast System (HAFS). The goal is to improve HAFS’s skill in operational forecast of TCs at the current resolution and higher resolution in future.

Results: Improved HAFS’s forecasting skill by the 3D-SA TKE scheme

We have developed an all-weather condition 3D-SA TKE scheme for UFS models and implemented it in the latest operational HAFSv2. The scheme has several unique features. First, it includes anisotropic TKE shear production tensor and 3D TKE turbulent transport and pressure diffusion. Second, it is a ‘wet’ scheme formulated with thermodynamic variables conserved for the moist adiabatic processes to allow for an appropriate consideration of the buoyancy induced by clouds. Third, it considers non-local vertical mixing induced by large turbulent eddies generated by the bottom-up forcing associated with surface heating and top-down forcing associated with cloud-top radiative cooling and entrainment instability of shallow clouds. Fourth, it partially unifies the vertical turbulent mixing parameterization with the FV3 horizontal divergence damping based on the predicted TKE. Finally, scale awareness is realized by blending turbulent eddy diffusivity, non-local mixing fluxes, and TKE dissipation length scale at the LES and LSMS extremes using partition functions empirically formulated based on the ratio of model resolution to the scale of large turbulent eddies19. Detailed description of the scheme is presented in Methods. The research-to-operation transition of the scheme for the operational HAFS is now underway at the Environmental Modeling Center (EMC) of NOAA. A version of HAFSv2 with the 3D-SA TKE scheme is available at https://github.com/hafs-community/HAFS/tree/feature/hafs.v2.0.1_3dtke.

The model tested in this study is HAFSv2.0.1, which is the baseline for the upcoming HAFSv2.1 targeting to be implemented for operational forecast in July 2025. HAFSv2.0.1 utilizes the new operational Real Time Ocean Forecast System version 2.5 (RTOFSv2.5). To evaluate the performance of the 3D-SA TKE scheme in HAFS’s real-time operational forecast setting, which currently uses a horizontal resolution of 2 km, EMC conducted two sets of parallel experiments using HAFSv2.0.1, named HFXR and HFXT hereafter. They are exactly the same except that HFXR uses HAFS’s default 1D Eddy-Diffusivity-Mass-Flux (EDMF) TKE scheme21, whereas HFXT activates the 3D-SA TKE scheme. Each set of experiments simulated 11 Atlantic basin hurricanes and tropical storms of the 2024 season. Table 1 summarizes the names of storms, storm intensity, date and time of initial and last cycle, and number of cycles being simulated for each storm. In total, 290 cycles have been simulated for each set of experiments.

Table 1 Information about the storms simulated by HFXR and HFXT
Full size table

Figure 2 compares the performance of storm track prediction by HFXR and HFXT. The spread of the forecast errors with respect to the Best-Track data is measured by the standard error of the mean (SEM) defined as \(\sqrt{\frac{{\sum }_{i=1}^{N}{({x}_{i}-\bar{x})}^{2}}{(N-1)\cdot N}}\), where \(x\) denotes a generic variable and \(N\) is the number of samples. The statistics show that there is a clear improvement in storm track forecast by the 3D-SA TKE scheme after 72 h in lead time. The decomposition of the total track error into the along- and cross-track components shows that the improvement mainly comes from the cross-track forecast. There are multiple factors that may lead to the along- and cross-track errors including steering flow, convection patterns, and model limitations that fail to capture the complex interaction between the ambient flow and storm convection22,23. The nearly identical along-track errors produced by HFXR and HFXT indicate that the two sets of experiments simulated similar storm translation speeds, which are mainly governed by the ambient steering flow. The reduction of cross-track errors in HFXT may be attributed to the impact of 3D turbulent transport processes on the distribution and intensity of eyewall and rainband convection.

Fig. 2: Statistics of track errors of all simulation cycles.
figure 2

Statistics of a total storm track errors, b absolute cross-track errors, and c absolute along-track errors with respect to the Best-Track data for all 290 simulated cycles by the two sets of experiments (HFXR and HFXT). The light transparent shades indicate the standard error of the mean (SEM). The numbers in cyan color at the bottom indicate the number of simulated cycles at different forecast lead times.

Full size image

Figure 3 compares the performance of storm intensity prediction by HFXR and HFXT. The absolute errors of both the maximum surface wind speed (VMAX) and the minimum storm central pressure (PMIN) are reduced by the 3D-SA TKE scheme in the late forecast lead times between 96 and 120 h. A marginal reduction in absolute errors of VMAX by the 3D-SA TKE scheme is also seen for lead times between 36 and 72 h. The most prominent feature is that the 3D-SA TKE scheme reduces the negative bias of VMAX of HFXR consistently for most of the forecast lead times (Fig. 3b).

Fig. 3: Statistics of intensity errors of all simulation cycles.
figure 3

a absolute errors of the predicted maximum surface wind speed (VMAX), b error biases with sign of VMAX, c absolute errors of the predicted minimum storm central pressure (PMIN) with respect to the Best-Track data. The light color shades denote SEM. Numbers in cyan color at the bottom indicate the number of simulated cycles at different forecast lead times.

Full size image

To better understand the improvement of track and intensity forecast by the 3D-SA TKE scheme, we classified the total simulation cycles into strong and weak categories based on storm’s initial intensity \(\ge\) 64 knots (i.e., hurricane strength) or < 64 knots (i.e., tropical storm or weaker systems). Of the total 290 simulation cycles, 141 and 149 cycles are classified into the strong and weak categories, respectively. The results show that the improvements in storm track and intensity forecast shown in Figs. 2 and 3 mainly come from the strong-category cycles. Figure 4 shows the statistics of forecast errors of track, VMAX, and PMIN for the strong category cycles, where improvements in track forecast (mainly the cross-track forecast) are seen after 72 h in lead time, and the improvements in intensity forecast (both absolute errors in VMAX and PMIN) are shown as early as 36 h in lead time. The paired Student’s T-tests show that the improvements in track forecasts, mainly the cross-track forecasts (Fig. 4a, b), by the 3D-SA TKE scheme are significant at 95% confidence level for lead time beyond 72 h, and the improvements in intensity forecasts measured by the error biases with sign of VMAX and PMIN (Fig. 4d, f) are significant at 95% confidence level for lead time beyond 84 h. Improvements in intensity forecast measured by the absolute errors of VMAX and PMIN appear to be less significant, but small P-values of ~0.05–0.1 are seen in lead time of ~36–60 and ~100–120 h for VMAX (Fig. 4c), and ~48–72 and ~96–120 h for PMIN (Fig. 4e), respectively.

Fig. 4: Statistics of forecast errors of simulation cycles with hurricane strength.
figure 4

a total track errors, b absolute cross- and along-track errors, c absolute errors of VMAX, d error biases of VMAX with sign, e absolute errors of PMIN, and f error biases of PMIN with sign for simulation cycles with intensity \(\ge\) \(64\) knots. Light color shades denote SEM. Green and gray lines scaled to the axes on the right indicate the P-values of Student’s t-tests of the forecast errors by HFXR and HFXT. Numbers in cyan color at the bottom indicate the number of simulated cycles at different forecast lead times.

Full size image

In contrast, no apparent improvement or downgrade of track and intensity forecast is seen for the weak category cycles (Supplementary Note 3 and Supplementary Fig. 3). The Student’s t-tests show that the differences of track and intensity forecasts by HFXR and HFXT are not statistically significant. The results shown in Fig. 4 and Supplementary Fig. 3 suggest that an appropriate account for 3D turbulent processes becomes increasingly important for forecasting stronger TCs. The reason is that large turbulent eddies generated by shear and buoyancy production become more energetic as a storm intensifies. Combining with the enhanced lateral contrast across the eyewall, rainbands, and moat in between, these energetic eddies result in stronger inter-connected 3D turbulent transport in the TC inner core.

We now examine the individual storms simulated by HFXR and HFXT. Figure 5 shows the composite storm intensities and tracks of Hurricanes Kirk (2024) and Leslie (2024) from multiple simulation cycles by HFXR and HFXT, respectively, compared with the Best-Track data. For Kirk (2024) and Leslie (2024), the 3D-SA TKE scheme reduces the intensity forecast error, in particular, HFXT captures the second intensity peak of Leslie (2024) to some extent after a temporary weakening, which is largely missed by HFXR. The 3D-SA TKE scheme also improves the storm track forecast of Kirk (2024) slightly as the TC vortex turns north and northeast in the later part of the simulations. In the case of Leslie (2024), the large spread in track forecast by HFXR in the late simulation hours (indicated by the large blue error bars) is substantially reduced in the HFXT simulations.

Fig. 5: HAFS’s prediction of hurricanes Kirk (2024) and Leslie (2024).
figure 5

a–f Composite storm intensity (VMAX and PMIN) and track of Hurricanes Kirk (2024) and Leslie (2024) from multiple simulation cycles by HFXR and HFXT, respectively, compared with the Best-Track data. The light color shades in the left and middle columns are SEM of the predicted intensity errors. The vertical and horizontal bars in the right column denote the standard deviations of the predicted storm center latitudes and longitudes with respect to the Best-Track data. The numbers in cyan color indicate the number of simulation cycles at different forecast lead times.

Full size image

To further examine the differences of simulated TC vortices between HFXR and HFXT, Fig. 6 shows the satellite images of 37 horizontally-polarized (37H) Special Sensor Microwave Imager/Sounder (SSMIS) brightness temperatures TB (K) of Hurricanes Kirk, Leslie, and Milton at an instant time during their intensification stage compared with the corresponding simulated 37H SSMIS TB (K) at the closest time to satellite images by HFXR and HFXT, respectively. Since the satellite observations have a low resolution of ~37\({\boldsymbol{\times }}\)28 km and operational HAFS has a resolution of 2 km, an apple-to-apple comparison is not possible. Nonetheless, both HAFS simulations qualitatively reproduced the observed magnitudes of TB and the global asymmetric structure of TC including the major spiral rainband features for all three storms. A noticeable difference between HFXR and HFXT is that the latter produced a somewhat more well-defined eyewall than the former and the structure in HFXT is closer to observations.

Fig. 6: Simulated TC structure compared with satellite observations.
figure 6

a, d, g Satellite images of 37 horizontally-polarized (37H) SSMIS brightness temperatures TB (K) of Hurricanes Kirk, Leslie, and Milton at an instant time during their intensification (curtesy to the Naval Research Laboratory Marine Meteorology Division). b, c, e, f, h, i The corresponding simulated 37H SSMIS TB (K) at the closest time to satellite images from HFXR and HFXT, respectively.

Full size image

Since there were reconnaissance flights flying into Milton (2024), the Milton simulations may be further evaluated against the airborne Tail Doppler Radar retrievals. Figure 7 compares the radar retrieved horizontal distribution of tangential wind at 1-km altitude and radius-height structure of azimuthal-mean radar reflectivity and tangential wind at 01Z Oct. 07, 2024, with the corresponding simulated fields from HFXR and HFXT. Note that the simulated tangential wind and reflectivity are the composites of 5 simulation cycles (2024100518–2024100618) at 00Z Oct. 07, 2024. Both HFXR and HFXT reasonably reproduced the observed magnitude of 1-km tangential wind and the size of TC vortex, but HFXT simulates a slightly better asymmetric structure of tangential winds. It is evident from the peak tangential wind in the west quadrant, which shifts slightly northward in HFXT consistent with the radar observations, but it shifts southward in HFXR. While both HAFS simulations appear to underestimate the observed vertical and radial extent of azimuthal-mean reflectivity and the vertical extent of azimuthal-mean tangential wind, HFXT produced a slightly better radius-height structure of azimuthal-mean reflectivity and tangential wind than HFXR in terms of vortex depth.

Fig. 7: Simulated Milton (2024) compared with Doppler radar observations.
figure 7

a, d, g Horizontal distribution of tangential wind at 1-km altitude and radius-height structure of azimuthal-mean reflectivity and tangential wind of Hurricane Milton (2024) at 01Z Oct. 7, 2024, from airborne tail Doppler radar retrievals. b, c, e, f, h, i The corresponding fields at 00Z Oct. 7, 2024, from HFXR and HFXT.

Full size image

Finally, it is worthwhile examining the predicted TKE from a scheme formulated with TKE. Figure 8 compares the azimuthal-mean predicted TKE averaged over the intensification period of Hurricane Kirk (00Z Oct. 2–00Z Oct. 4, 2024) from the simulation cycle of 2024092912 by HFXR and HFXT. Large TKE occurs in the low troposphere in the vicinity of the eyewall. The magnitude of simulated TKE is comparable to that retrieved from airborne radars12,13. Compared with HFXR, HFXT produces the larger TKE near the eyewall in the low troposphere where the horizontal shear is the most prominent. The large difference of TKE between HFXR and HFXT occurs near the boundary layer top (Fig. 8d), which is consistent with the LES result (Supplementary Fig. 2), suggesting that the 3D-SA TKE scheme captures the 3D turbulent transport processes in the vicinity of the eyewall. We also examined the sensitivity of the predicted TKE to model resolution in idealized HAFS simulations. The results show that the 3D-SA TKE scheme captures the eyewall TKE in the upper troposphere resembling that of Doppler radar retrievals12,13 as the model grid spacing reduces to 1 km (Supplementary Note 4 and Supplementary Fig. 4).

Fig. 8: Simulated TKE from Hurricane Kirk (2024).
figure 8

a, b Radius-height structure of azimuthal-mean simulated TKE averaged over 00Z Oct. 2–00Z Oct. 4, 2024, from the simulation cycle 20240912 of Hurricane Kirk (2024) by HFXR and HFXT, respectively. c Difference of azimuthal-mean TKE between HFXT and HFXR. d Vertical profiles of azimuthal-mean TKE at the radius of 50 km from HFXR and HFXT.

Full size image

As stated previously, the current operational HAFS uses a horizontal grid resolution of 2 km. To examine the sensitivity of the 3D-SA TKE scheme to model horizontal resolution, using HAFSv2.0.1 in the real-time operational forecasting setting except for horizontal resolution, we simulated Hurricane Milton (2024) with different horizontal resolutions of 800 m, 1 km, 1.5 km, 2 km, 4 km, and 8 km. For each resolution, 21 simulation cycles from 2024100512 to 2024101012 with an interval of 6 hours were conducted. Due to the limited computational resources, the 800-m resolution simulations can only be run for ~54 – 60 forecast hours. Therefore, comparisons among the sensitivity experiments with different resolutions were only made up to 54 hours, although low resolution simulations were run for 126 hours. Figure 9 compares the errors of the simulated storm track and intensity with respect to the Best-Track data. The statistics show that the higher resolution simulations (800 m, 1 km, and 1.5 km) overall produce smaller track and intensity errors than the lower resolution simulations (2 km, 4 km, and 8 km) although with exceptions. For example, the 8-km resolution simulations produce the smallest cross-track error, but they also produce the largest along-track error, which dominates the cross-track error.

Fig. 9: Sensitivity to model horizontal resolution.
figure 9

Statistics of a total track errors, b absolute along-track errors, c absolute cross-track errors, d absolute errors of VMAX, e error biases of VMAX with sign, and f absolute errors of PMIN as the function of forecast lead time for simulations of Hurricane Milton (2024) at different horizontal resolutions of 800 m, 1 km, 1.5 km, 2 km, 4 km, and 8 km. Numbers in cyan color at the bottom indicate the number of simulated cycles at different forecast lead times.

Full size image

Two things need to be acknowledged here. First, the results could be affected by the small sample size from simulations of one storm. Therefore, there is an uncertainty of the sensitivity shown in Fig. 9. Second, evaluation of the sensitivity of a turbulence scheme to model resolution under complex TC conditions is a complicated problem. Unlike the dry CBL whose characteristics are largely determined by the turbulent mixing processes, the evolution of a TC and its inner-core structure also depend strongly on cloud microphysics and other processes in addition to turbulent mixing. As model resolution increases, fine scale cloud processes resolved by the model and possibly improved representation of microphysics at the finer grids can all affect TC simulations. How to appropriately separate the effect of turbulent mixing from other complications resulting from moist convection is an important issue that needs to be addressed. We will tackle the problem in our future study. Nonetheless, the sensitivity results shown in Fig. 9 at least suggest that the 3D-SA TKE scheme functions well in TC simulations for a range of horizontal resolutions.

Discussion

How to realistically parameterize SGS anisotropic 3D turbulent mixing in the gray zone is a longstanding problem of NWP. Previous attempts for developing 3D turbulence parameterization in the gray zone mainly focused on the dry CBL in horizontally homogenous conditions19,24. In this study, we have developed a 3D-SA TKE scheme applicable to all conditions including weather extremes, such as TCs, for UFS operational models with the FV3 dynamic core and implemented it in HAFS. The goal is to enhance operational HAFS’s skill in predicting storm track and intensity of TCs at the current resolution and higher resolutions in future. The developed 3D-SA TKE scheme has been tested in the latest version of HAFS (HAFSv2.0.1) targeting operational implementation of the 2025 season. EMC conducted two sets of parallel experiments that simulated 11 Atlantic basin storms of the 2024 season in HAFS real-time operational forecast setting totaling 290 cycles. These two experiments are exactly the same except that one uses HAFS’s default 1D EDMF TKE scheme, and the other activates the 3D-SA TKE scheme. Comparison of the simulated storm tracks and intensities with the Best-Track data shows that the 3D-SA TKE scheme improves HAFS forecast skill by reducing the forecasting errors of cross-track, VMAX, and PMIN in day 4 and day 5 of forecast lead time. Classifying the total simulation cycles into the strong and weak categories based on the intensity threshold of 64 knots (i.e., hurricane strength) reveals that the improvements in forecasting skill mainly come from the strong category where the improvements are statistically significant at 95% confidence level according to the Student’s t-tests, suggesting that an appropriate account for the 3D turbulent transport processes becomes increasingly important for track and intensity forecasts as a storm intensifies. Comparisons with the satellite and airborne radar observations also show that the HAFS with the 3D-SA TKE scheme produces better vortex inner core structures during the storm intensification than the default HAFS. In short, validations of model’s performance in HAFS operational setting show promising results of the 3D-SA TKE scheme in enhancing operational HAFS’s skill in storm track and intensity prediction. The transition from research to operation is now underway at EMC of NOAA. Finally, we note that although the 3D-SA TKE scheme has been developed within HAFS, its implementation in the FV3 dynamic core and the Common Community Physics Package (CCPP25) makes it ready to be used in other UFS models. In addition to the research-to-operation transition of the scheme in HAFS, EMC is now also testing the 3D-SA TKE scheme in operational Global Forecast System (GFS) model.

For future work, we will focus on the following areas. First, the main task of this study was to develop a 3D-SA TKE scheme for UFS models and examine its performance in HAFS operational setting that uses a 2-km resolution for the nested ___domain. At such a resolution, the benefit of using a 3D-SA turbulence scheme may be limited. It is expected that the 3D turbulence parameterization will have a stronger impact on TC simulations as model resolution increases. Although in this study we conducted experiments to examine the sensitivity of the 3D-SA TKE scheme to model resolution, only one storm has been tested with limited cycles. We are planning to continue exploring the sensitivity of HAFS simulations to model grid resolution when more computational resources are available. Second, as stated in Methods, HAFS uses different turbulence closures to parameterize SGS turbulent mixing: the classic K-closure for handling vertical mixing closed by eddy diffusivity and divergence damping for treating horizontal mixing closed by damping coefficients. In this study, we integrated the second order divergence damping with the vertical turbulent mixing parameterization based on the predicted TKE. We will continue working with scientists at the Geophysical Fluid Dynamics Laboratory (GFDL) of NOAA to explore dynamically consistent methods to apply the 3D-SA TKE-based damping to physical quantities beyond the horizontal winds within FV3. Third, recent studies showed that the improvements in 1D PBL scheme enhance the skill of HAFS-B in predicting TCs26,27. This motivates us to further examine the performance of the 3D-SA TKE scheme in HAFS-B simulations in our future study. Finally, current numerical models consist of a dynamic core and a separate physics package constructed in 1D vertical column framework. Such a model infrastructure works well for coarse-resolution models. However, as model resolution increases, the 3D physical processes become increasingly important, which raises a concern if the coupling method by transferring data between the separated dynamic core and 1D physics modules can realistically represent the dynamics-physics interaction in the real atmosphere. In the case of our 3D-SA TKE scheme, it involves both FV3 dynamic core and 1D CCPP. The data transfer between them is a potential obstacle hindering model efficiency, accurate implementation of physical processes, and tight coupling between dynamics and physics. One method to solve the problem is to integrate physical processes directly within the dynamic core. One notable success has been with the GFDL microphysics28,29, which is implemented directly within the dynamic core30, allowing for a tight integration between dynamic and thermodynamic processes. As indicated by Eq. (1) in Methods, the tendencies induced by turbulence are directly involved in the momentum budget, thus, it requires a tight coupling between the resolved dynamic processes and SGS turbulence processes. This issue may be addressed by integrating the 3D turbulent mixing parameterization directly within the FV3 dynamic core. We will explore the possibility in our future studies.

Methods: A 3D-SA TKE scheme applicable to the gray zone

The grid-box mean governing equations of momentum may be expressed as,

$$\frac{\partial {\bar{u}}_{i}}{\partial t}=-{\bar{u}}_{j}\frac{\partial {\bar{u}}_{i}}{\partial {x}_{j}}-\frac{\partial \bar{p}}{\bar{\rho }\partial {x}_{i}}-2{\epsilon }_{{ijk}}{\varOmega \eta }_{j}{\bar{u}}_{k}+{\delta }_{i3}g-\frac{\partial \overline{{u}_{i}^{{\prime} }{u}_{j}^{{\prime} }}}{\partial {x}_{j}},$$
(1)

where \(u\) is the velocity component, subscript indices \(i,j,{k}=1,\,2,\,3\) indicate \(x,{y},{z}\) directions in a local Cartesian coordinate, \(p\) is the pressure, \(\rho\) is the air density, \(g\) is the gravity, \(\varOmega\) is the Earth’s angular velocity, \({\eta }_{j}=[0,\,\cos \phi ,\,\sin \phi\)], \(\phi\) is the latitude, \({\delta }_{{ij}}\) is the Kronecker delta function, \({\epsilon }_{{ijk}}\) is the Levi-Civita function, overbar and prime indicates the grid-box mean and perturbation away from the mean. \(-\frac{\partial \overline{{u}_{i}^{{\prime} }{u}_{j}^{{\prime} }}}{\partial {x}_{j}}\) represents the tendencies that must be determined parametrically by a turbulence scheme. At the LES resolution, only isotropic small eddies need to be parameterized, therefore, turbulent flux is often assumed to follow the local down-gradient transport theory3,

$$-\overline{{u}_{i}^{{\prime} }{u}_{j}^{{\prime} }}={k}_{m({les})}{\bar{S}}_{{ij}},\,{\bar{S}}_{{ij}}=\left(\frac{\partial {\bar{u}}_{i}}{\partial {x}_{j}}+\frac{\partial {\bar{u}}_{j}}{\partial {x}_{i}}\right),\,i,j=1,\,2,\,3,$$
(2)

where \({k}_{m({les})}\) is the eddy diffusivity at the LES extreme, which is the same in all directions under the isotropic assumption. In the 3D Smagorinsky–Lilly model, \({k}_{m({les})}\) itself is assumed to be the function of local 3D gradients,

$${k}_{m({les})}={c}_{s}^{2}{l}_{({les})}^{2}\left|{\bar{S}}_{{ij}}\right|,$$
(3)

where \({c}_{s}\) is Smagorinsky constant. Under the assumption that the rate of energy cascade equals the rate of energy dissipation, it can be shown \({c}_{s}={(\frac{3{\pi }^{4/3}}{2}\alpha )}^{-3/4}\), where α is Kolmogorov constant3. In LES applications, \({c}_{s}\) often takes a value of ~0.1 – 0.28. \({l}_{({les})}\) is the turbulent mixing length at the order of model grid spacing.

Alternatively, \({k}_{m({les})}\) may be parameterized in terms of TKE, \(e=\tfrac{1}{2}\overline{{u}_{i}^{{\prime} }{u}_{i}^{{\prime} }}\), as,

$${k}_{m({les})}={c}_{m}{l}_{({les})}\sqrt{e},$$
(4)

where \({c}_{m}\) is a coefficient. The first 3D TKE scheme was developed to simulate the mixing layer topped by stratocumulus4. It is now known as the 3D Deardorff TKE model in which TKE is prognostically determined by the TKE budget equation,

$$\frac{\partial e}{\partial t}=-{\bar{u}}_{j}\frac{\partial e}{\partial {x}_{j}}-\overline{{u}_{\text{i}}^{{\prime} }{u}_{\text{j}}^{{\prime} }}\frac{\partial {\bar{u}}_{i}}{\partial {x}_{j}}+{\delta }_{i3}\frac{g}{\overline{{\theta }_{v}}}\overline{{{u}_{\text{i}}^{{\prime} }\theta }_{v}^{{\prime} }}-\frac{\partial }{\partial {x}_{i}\,}\left[\,\overline{{u}_{\text{i}}^{{\prime} }\left(e+\frac{{p}^{{\prime} }}{\bar{\rho }}\right)}\,\right]-\varepsilon ,$$
(5)

where \({\theta }_{v}\) is the virtual potential temperature and \(\varepsilon\) is the TKE dissipation rate parameterized as,

$$\varepsilon =\frac{{e}^{3/2}}{{l}_{({les})}/{c}_{\varepsilon }},$$
(6)

and dissipation coefficient \({c}_{\varepsilon }\) can be shown to be \({c}_{\varepsilon }=\pi {(\frac{2}{3\alpha })}^{\frac{3}{2}}\) 3.\(\frac{{l}_{({les})}}{{c}_{\varepsilon }}\) is known as the TKE dissipation length scale at the LES extreme\(.\)

In fact, Eq. (4) may be derived from the simplified TKE budget equation. Assuming that there is a balance between shear production and dissipation of TKE, Eq. (5) reduces to \(-\overline{{u}_{i}^{{\prime} }{u}_{j}^{{\prime} }}\frac{\partial {\bar{u}}_{i}}{\partial {x}_{j}}=\varepsilon\). Inserting Eq. (2), Eq. (3), and Eq. (6) into this simplified equation, it yields Eq. (4) with \({c}_{m}={c}_{s}^{4/3}{c}_{\varepsilon }^{1/3}=\frac{1}{\pi }{(\frac{2}{3\alpha })}^{3/2}\). One limitation of the 3D Smagorinsky–Lilly model is that the effect of static stability on turbulent mixing is not considered. To solve the problem, stability correction is often added to the 3D Smagorinsky-Lilly model to account for the complication from the thermal stratification31. The 3D Deardorff TKE model, on the other hand, successfully included the stability effects on SGS turbulent mixing since the TKE budget includes buoyancy production or suppression of TKE. Moreover, a TKE scheme does not need PBL height information and can be applied to both deep and shallow convection regimes provided that TKE production, suppression, transport, and dissipation can be appropriately determined. The eddy diffusivity for heat and moisture may be taken as4,

$$\left\{\begin{array}{ll}{k}_{h\left({les}\right)}={k}_{q\left({les}\right)}=\left(1+2\frac{{l}_{({les})}}{\Delta s}\right){k}_{m\left({les}\right)},\,\Delta s={(\Delta x\Delta y\Delta z)}^{1/3}\,\\\;\; {l}_{({les})}=0.76\sqrt{e}{(\frac{g}{\overline{{\theta }_{l}}}\frac{\partial \overline{{\theta }_{l}}}{\partial z})}^{-1/2},\,\frac{\partial \overline{{\theta }_{l}}}{\partial z} > 0\,\\ \;\;{l}_{({les})}=\Delta s,\,\frac{\partial \overline{{\theta }_{l}}}{\partial z}\le 0\,\end{array}\right.,$$
(7)

where \(\bar{{\theta }_{l}}\) is the liquid water potential temperature.

Applying Eq. (2), the 3D shear production of TKE in Eq. (5) may be expressed in terms of eddy diffusivity and 3D wind shear as,

$$-\overline{{u}_{\text{i}}^{{\prime} }{u}_{\text{j}}^{{\prime} }}\frac{\partial {\bar{u}}_{i}}{\partial {x}_{j}}={k}_{m\left({les}\right)}\left\{2\left[{(\frac{\partial \bar{u}}{\partial x})}^{2}+{(\frac{\partial \bar{v}}{\partial y})}^{2}+{(\frac{\partial \bar{w}}{\partial z})}^{2}\right]+{\left(\frac{\partial \bar{u}}{\partial y}+\frac{\partial \bar{v}}{\partial x}\right)}^{2}+{\left(\frac{\partial \bar{u}}{\partial z}+\frac{\partial \bar{w}}{\partial x}\right)}^{2}+{\left(\frac{\partial \bar{v}}{\partial z}+\frac{\partial \bar{w}}{\partial y}\right)}^{2}\right\}$$
(8)

In a 1D TKE scheme, only vertical shear \(\left\{{(\frac{\partial \bar{u}}{\partial z})}^{2}+{(\frac{\partial \bar{v}}{\partial z})}^{2}\right\}\) is considered. However, Eq. (8) cannot be extended to the gray zone since eddy diffusivity is no longer isotropic in the gray zone.

To appropriately account for the directionally different 3D turbulent mixing induced by the anisotropic eddies in the gray zone, the fundamental turbulence closure needs to be modified. To allow for anisotropic turbulent mixing, Eq. (2) may be generalized as,

$$-\overline{{u}_{i}^{{\prime} }{u}_{j}^{{\prime} }}={k}_{m}\left(i,j\right)\frac{\partial {\bar{u}}_{i}}{\partial {x}_{j}}+{k}_{m}\left(j,i\right)\frac{\partial {\bar{u}}_{j}}{\partial {x}_{i}},$$
(9)

where \({k}_{m}\left(i,j\right)\) and \({k}_{m}\left(j,i\right)\) are the shear-direction dependent eddy diffusivity. Note that \({k}_{m}\left(i,j\right)\) does not necessarily equal \({k}_{m}\left(j,i\right)\) due to the anisotropic nature of eddies. For example, for \(i=1,j=3\), Eq. (9) becomes \(-\overline{{u}^{{\prime} }{w}^{{\prime} }}={k}_{m}\left(x,z\right)\frac{\partial \bar{u}}{\partial z}+{k}_{m}\left(z,x\right)\frac{\partial \bar{w}}{\partial x}\), where \({k}_{m}\left(x,z\right)\) denotes the vertical eddy diffusivity associated with the vertical gradient of x-direction wind component, whereas \({k}_{m}\left(z,x\right)\) is the horizontal eddy diffusivity associated with the x-direction gradient of vertical velocity. For isotropic conditions, \({k}_{m}\left(i,j\right)={k}_{m}\left(j,i\right)={k}_{m({les})}\) and Eq. (9) reduces to Eq. (2). Assuming that all horizontal eddy diffusivities associated with horizontal gradients are the same, and all vertical eddy diffusivities associated with vertical gradients are the same, i.e.,

$$\left\{\begin{array}{ll}{k}_{m}\left(x,x\right)={k}_{m}\left(y,x\right)={k}_{m}\left(z,x\right)={k}_{m}\left(x,y\right)={k}_{m}\left(y,y\right)={k}_{m}\left(z,y\right)={k}_{m}^{(h)}\\ {k}_{m}\left(x,z\right)={k}_{m}\left(y,z\right)={k}_{m}\left(z,z\right)={k}_{m}^{(v)}\,\end{array}\right.,$$
(10)

where \({k}_{m}^{(h)}\) and \({k}_{m}^{(v)}\) are the horizontal and vertical eddy diffusivity in the gray zone, respectively. Then, using Eq. (9) and Eq. (10), the 3D shear production becomes,

$$\begin{array}{c}-\overline{{u}_{i}^{{\prime} }{u}_{j}^{{\prime} }}\frac{\partial {\bar{u}}_{i}}{\partial {x}_{j}}={k}_{m}^{(h)}\left\{2\left[{(\frac{\partial \bar{u}}{\partial x})}^{2}+{(\frac{\partial \bar{v}}{\partial y})}^{2}\right]+{\left(\frac{\partial \bar{u}}{\partial y}+\frac{\partial \bar{v}}{\partial x}\right)}^{2}+\left[{(\frac{\partial \bar{w}}{\partial x})}^{2}+{(\frac{\partial \bar{w}}{\partial y})}^{2}\right]+\left(\frac{\partial \bar{w}}{\partial x}\frac{\partial \bar{u}}{\partial z}+\frac{\partial \bar{w}}{\partial y}\frac{\partial \bar{v}}{\partial z}\right)\right\}\\ +{k}_{m}^{(v)}\left\{2{(\frac{\partial \bar{w}}{\partial z})}^{2}+\left[{(\frac{\partial \bar{u}}{\partial z})}^{2}+{(\frac{\partial \bar{v}}{\partial z})}^{2}\right]+\left(\frac{\partial \bar{w}}{\partial x}\frac{\partial \bar{u}}{\partial z}+\frac{\partial \bar{w}}{\partial y}\frac{\partial \bar{v}}{\partial z}\right)\right\}.\end{array}$$
(11)

For isotropic turbulence, \({k}_{m}^{(h)}\) and \({k}_{m}^{(v)}\) becomes \({k}_{m\left({les}\right)}\), and Eq. (11) reduces to Eq. (8). \({k}_{m}^{(h)}\) and \({k}_{m}^{(v)}\) in the gray zone may be obtained by blending the eddy diffusivities at the LES and LSMS extremes19,

$$\left\{\begin{array}{c}{k}_{m}^{(h)}=\left(1-{f}_{p}\right){k}_{m({les})}+{f}_{p}{k}_{m({lsms})}^{(h)}\\ {k}_{m}^{(v)}=\left(1-{f}_{p}\right){k}_{m({les})}+{f}_{p}{k}_{m({lsms})}^{(v)}\end{array}\right.,$$
(12)

where \({k}_{m({lsms})}^{(h)}\) and \({k}_{m({lsms})}^{(v)}\) are the horizontal and vertical eddy diffusivities at the LSMS extreme. Since the 3D-SA TKE scheme will asymptotically approach to a 1D PBL scheme at the LSMS extreme, and in our case, it is designed for HAFS, \({k}_{m({lsms})}^{(v)}\) is taken as that of HAFS default 1D EDMF TKE scheme21. However, since the FV3 dynamic core does not use classic K-closure but employs the so-called divergence damping to parameterize horizontal momentum diffusion, \({k}_{m({lsms})}^{(h)}\) is not available in HAFS (Details of how to incorporate the 3D-SA TKE scheme with FV3 are provided in Section 4.5). For this reason, \({k}_{m({lsms})}^{(h)}\) is simply taken as,

$${k}_{m({lsms})}^{(h)}={c}_{{mh}}{A}^{1/2}\sqrt{e},$$
(13)

where \({c}_{{mh}}\) is an adjustable empirical coefficient, here taken as 0.1. \(A\) is the horizontal area of a grid box. One reason to choose Eq. (13) is that both horizontal and vertical diffusion can be parameterized consistently based on the predicted TKE. The horizontal eddy diffusivity for heat and moisture at the LSMS extreme, then, is taken as,

$${k}_{h({lsms})}^{(h)}={k}_{q({lsms})}^{(h)}={k}_{m({lsms})}^{(h)}{P}_{r}^{-1},$$
(14)

where \({P}_{r}\) is the turbulent Prandtl number.

The SA feature of the scheme is realized though a partition function \({f}_{p}={f}_{p}(\frac{\sqrt{A}}{{S}_{{le}}})\), where \({S}_{{le}}\) denotes the scale of large eddies. When the model grid spacing \(\sqrt{A}\) is much smaller than \({S}_{{le}}\) (LES extreme), the unresolved small eddies become isotropic, \({f}_{p}=0\), and the directionally different eddy diffusivities \({k}_{m}^{(h)}\) and \({k}_{m}^{(v)}\) become isotropic eddy diffusivity \({k}_{m({les})}\). Conversely, when the model grid spacing \(\sqrt{A}\) is much larger than \({S}_{{le}}\), \({f}_{p}=1\) and \({k}_{m}^{(h)}\) and \({k}_{m}^{(v)}\) reduce to \({k}_{m({lsms})}^{(h)}\) and \({k}_{m({lsms})}^{(v)}\) at the LSMS extreme. For our 3D-SA TKE scheme, \({S}_{{le}}\) is set to the larger value between the diagnosed PBL height and the height where TKE decreases to 50% of the maximum TKE in a vertical column. This is based on the consideration that large eddies in the boundary layer are constrained by the capping inversion if it exists. However, when an inversion does not exist like the deep convective eyewall and rainbands (Fig. 1), the PBL height diagnosed from the bulk Richardson number21 may underestimate the scale of large eddies, therefore, we use TKE as an additional criterion for estimating the scales of large eddies. The criterion of 50% of the maximum TKE is chosen based on the sensitivity tests (Supplementary Note 5 and Supplementary Fig. 5). For very stable conditions, \({S}_{{le}}\) takes the mixing length computed in the 1D TKE scheme. The empirical formula of \({f}_{p}(\frac{\sqrt{A}}{{S}_{{le}}})\) derived from LESs19 is adopted in our scheme. The same blending as Eq. (12) is, then, applied to the eddy diffusivity of heat and moisture to obtain \({k}_{h}^{(h)}\), \({k}_{q}^{(h)}\), \({k}_{h}^{(v)}\), and \({k}_{q}^{(v)}\) in the gray zone.

The eddy transport of TKE and pressure diffusion, \(-\frac{\partial }{\partial {x}_{i}\,}\left[\,\overline{{u}_{i}^{{\prime} }\left(e+\frac{{p}^{{\prime} }}{\bar{\rho }}\right)}\,\right]\), is another term in the TKE budget that needs to be treated anisotropically in the gray zone. Following the anisotropic down-gradient transport of Eq. (9), the term may be parameterized as,

$$-\frac{\partial }{\partial {x}_{i}\,}\left[\,\overline{{u}_{\text{i}}^{{\prime} }\left(e+\frac{{p}^{{\prime} }}{\bar{\rho }}\right)}\,\right]\,=\frac{\partial }{\partial x}\left(2{k}_{q}^{\left(h\right)}\frac{\partial e}{\partial x}\right)+\frac{\partial }{\partial y}\left(2{k}_{q}^{\left(h\right)}\frac{\partial e}{\partial y}\right)+\frac{\partial }{\partial z}\left(2{k}_{q}^{\left(v\right)}\frac{\partial e}{\partial z}\right).$$
(15)

For isotropic small eddies, \({k}_{q}^{\left(h\right)}\) and \({k}_{q}^{\left(v\right)}\) become \({k}_{q({les})}\) at the LES extreme and Eq. (15) reduces to that of the 3D Deardorff TKE model4. TKE is treated as a tracer in HAFS similar to moisture and hydrometeors, thus, \({k}_{q}^{\left(h\right)}\) and \({k}_{q}^{\left(v\right)}\) are used in Eq. (15). With the known bottom and top boundary conditions of TKE, \(\frac{\partial }{\partial z}\left(2{k}_{q}^{\left(v\right)}\frac{\partial e}{\partial z}\right)\) is solved implicitly via a tridiagonal solver. However, since the lateral boundary conditions of TKE are unknown due to the parallel calculation depending on the assigned computer nodes, \(\scriptstyle\frac{\partial }{\partial x}\left(2{k}_{q}^{\left(h\right)}\frac{\partial e}{\partial x}\right)+\frac{\partial }{\partial y}\left(2{k}_{q}^{\left(h\right)}\frac{\partial e}{\partial y}\right)\) is calculated as an explicit tendency term and added into the TKE budget equation.

Finally, the parameterization of TKE dissipation in the gray zone takes the form of32,

$$\varepsilon ={e}^{3/2}/{l}_{\varepsilon },$$
(16)

where \({l}_{\varepsilon }\) is the TKE dissipation length scale in the gray zone, which may be parameterized by blending the corresponding length scales at the LES and LSMS extremes as,

$${l}_{\varepsilon }=\left(1-{f}_{\varepsilon p}\right)\frac{{l}_{\left({les}\right)}}{{c}_{\varepsilon }}+{f}_{\varepsilon p}\frac{{l}_{({lsms})}^{(v)}}{{c}_{\varepsilon ({lsms})}},\,{f}_{\varepsilon p}={f}_{\varepsilon p}\left(\frac{\sqrt{A}}{{S}_{{le}}}\right),$$
(17)

where \({f}_{\varepsilon p}\) is the partition function for dissipation coefficients19. For \({f}_{\varepsilon p}=0\,{\rm{or}}\,1\), Eq. (17) reduces to \(\frac{{l}_{\left({les}\right)}}{{c}_{\varepsilon }}\) and \(\frac{{l}_{({lsms})}^{(v)}}{{c}_{\varepsilon ({lsms})}}\) at the LES and LSMS extremes, respectively. Since \({c}_{\varepsilon }\) considers isotropic 3D TKE dissipation and \({c}_{\varepsilon ({lsms})}\) considers 1D TKE dissipation, Eq. (17) provides a way to account for anisotropic 3D TKE dissipation in the gray zone as \(\frac{\sqrt{A}}{{S}_{{le}}}\) changes. The vertical mixing length \({l}_{({lsms})}^{(v)}\) and TKE dissipation coefficient \({c}_{\varepsilon ({lsms})}\) take those of 1D EDMF TKE scheme21 since the 3D-SA scheme asymptotically approaches a 1D PBL scheme at the LSMS extreme.

The 3D-SA TKE scheme is designed for UFS operational models applicable to all weather conditions. Latent heating is an important mechanism for the buoyancy production of TKE in the clouds. Studies show that the turbulent transport induced by cloud processes above the boundary layer in the eyewall and rainbands plays an important role in TC intensification13,33. The effect of cloud induced buoyancy on TKE generation may be included in the scheme using the Brunt-Väisälä frequency \({N}^{2}\) for the saturated atmosphere containing liquid- and solid-phase hydrometeors. For a “dry” scheme formulated by the thermodynamic variables not conserved for moist adiabatic process, such as, temperature, potential temperature, and water vapor mixing ratio, a stability correction is needed by explicitly including the effects of liquid- and solid-phase hydrometeors in the formula of \({N}^{2}\)33. Alternatively, \({N}^{2}\) can be calculated using the thermodynamic variables conserved for moist adiabatic process, such as liquid-ice static energy \({H}_{{li}}={c}_{{pd}}T+{gz}-L\left({q}_{l}+{q}_{i}\right)-{L}_{f}{q}_{i}\) and total water mixing ratio \({q}_{t}\),

$${N}^{2}={b}_{h}\frac{\partial {\bar{H}}_{{li}}}{\partial z}+{b}_{q}\frac{\partial {\bar{q}}_{t}}{\partial z},$$
(18)

where \({b}_{h}\) and \({b}_{q}\) are coefficients resulting from the derivations of moist thermodynamics21,34. Once \({N}^{2}\) is determined, the buoyancy production of TKE can be calculated as,

$$\frac{g}{\overline{{\theta }_{v}}}\overline{{{w}^{{\prime} }\theta }_{v}^{{\prime} }}=-{k}_{h}^{\left(v\right)}\frac{g}{\bar{{\theta }_{v}}}\frac{\partial \bar{{\theta }_{v}}}{\partial z}=-{k}_{h}^{\left(v\right)}{N}^{2}.$$
(19)

The non-local mixing at the LES extreme is typically negligible as the local down-gradient transport induced by small eddies dominates, but it becomes increasingly important as model grid resolution decreases into the gray zone. For the current 3D-SA TKE scheme, only vertical non-local mixing is considered. We are working on the methods to parameterize horizontal non-local mixing, which will be included in the future version of the scheme. The total vertical turbulent fluxes in the gray zone may be decomposed into local and non-local components as,

$$\left\{\begin{array}{c}\overline{{u}_{i}^{{\prime} }{w}^{{\prime} }}=-{k}_{m}^{(v)}\frac{\partial {\bar{u}}_{i}}{\partial z}+{\overline{{u}_{i}^{{\prime} }{w}^{{\prime} }}}^{{NL}}\\ \overline{{\phi }^{{\prime} }{w}^{{\prime} }}=-{k}_{h}^{(v)}\frac{\partial \bar{\phi }}{\partial z}+{\overline{{\phi }^{{\prime} }{w}^{{\prime} }}}^{{NL}}\end{array}\right.,$$
(20)

where \(\phi\) is a generic scalar. \(-{k}_{m}^{(v)}\frac{\partial {\bar{u}}_{i}}{\partial z}\) and \(-{k}_{h}^{(v)}\frac{\partial \bar{\phi }}{\partial z}\) are handled by the local down-gradient parameterization described previously. \({\overline{{u}_{i}^{{\prime} }{w}^{{\prime} }}}^{{NL}}\) and \({\overline{{\phi }^{{\prime} }{w}^{{\prime} }}}^{{NL}}\) represent the non-local vertical fluxes in the gray zone, which can be calculated by blending their counterparts at the LES and LSMS extremes19,

$$\left\{\begin{array}{ll}{\overline{{u}_{i}^{{\prime} }{w}^{{\prime} }}}^{{NL}}={f}_{{pnl}}\left(\frac{\sqrt{A}}{{S}_{{le}}}\right){\overline{{u}_{i}^{{\prime} }{w}^{{\prime} }}}_{({lsms})}^{{NL}},\,i=1,\,2\\ {\overline{{\phi }^{{\prime} }{w}^{{\prime} }}}^{{NL}}={f}_{{pnl}}\left(\frac{\sqrt{A}}{{S}_{{le}}}\right){\bar{{\phi }^{{\prime} }{w}^{{\prime} }}}_{({lsms})}^{{NL}}\,\end{array}\right.,$$
(21)

where \({f}_{{pnl}}\left(\frac{\sqrt{A}}{{S}_{{le}}}\right)\) is a partition function for non-local fluxes with \({f}_{{pnl}}=0\) and 1 at the LES and LSMS extremes respectively. Note that the non-local transport is assumed to be zero at the LES extreme. The non-local vertical fluxes at the LSMS extreme are determined using the mass-flux approach21,

$$\left\{\begin{array}{c}{\overline{{u}_{i}^{{\prime} }{w}^{{\prime} }}}_{({lsms})}^{{NL}}={M}_{{up}}\left({u}_{{iup}}-\bar{{u}_{i}}\right)+{M}_{{dw}}\left({u}_{{idw}}-\bar{{u}_{i}}\right)\,\\ {\overline{{\phi }^{{\prime} }{w}^{{\prime} }}}_{({lsms})}^{{NL}}={M}_{{up}}\left({\phi }_{{up}}-\bar{\phi }\right)+{M}_{{dw}}\left({\phi }_{{dw}}-\bar{\phi }\right)\,\end{array}\right.,$$
(22)

where \({M}_{{up}}\) and \({M}_{{dw}}\) are the convective mass flux associated with the coherent updrafts and downdrafts, and \({u}_{{iup}}\), \({u}_{{idw}}\), \({\phi }_{{up}}\), and \({\phi }_{{dw}}\) are the horizontal wind components and scalar associated with updrafts and downdrafts, respectively. The terms on the right-hand side of Eq. (22) represent the non-local fluxes induced by the coherent updrafts and downdrafts resulting from the bottom-up forcing associated with the surface heating and top-down forcing associated with the cloud-top radiative cooling and cloud-top entrainment instability of shallow stratocumulus, which are calculated following 1D EDMF TKE scheme21.

As stated previously, turbulence parameterization is to parametrically determine the tendencies induced by SGS turbulence in the grid-box mean governing equations. Using the momentum budget as an example, it is \(-\frac{\partial \overline{{u}_{i}^{{\prime} }{u}_{j}^{{\prime} }}}{\partial {x}_{j}}\). HAFS parameterizes the horizontal component of the tendency (\(-\frac{\partial \overline{{u}^{{\prime} }{u}^{{\prime} }}}{\partial x}-\frac{\partial \overline{{u}^{{\prime} }{v}^{{\prime} }}}{\partial y}\) and \(-\frac{\partial \overline{{v}^{{\prime} }{u}^{{\prime} }}}{\partial x}-\frac{\partial \overline{{v}^{{\prime} }{v}^{{\prime} }}}{\partial y}\)) in the FV3 dynamic core and the vertical component (\(-\frac{\partial \overline{{u}^{{\prime} }{w}^{{\prime} }}}{\partial z}\) and \(-\frac{\partial \overline{{v}^{{\prime} }{w}^{{\prime} }}}{\partial z}\)) in the 1D EDMF TKE scheme. In FV3, the so-called divergence damping is used to parameterize the horizontal momentum diffusion. The governing equations of horizontal winds are expressed as20,

$$\left\{\begin{array}{c}{\bar{u}}^{n+1}={\bar{u}}^{n}+\cdots +\frac{\delta }{\delta x}({\nu }_{d}{\nabla }^{2N}D)+\frac{\delta }{\delta x}({\nu }_{2d}D)\\ {\bar{v}}^{n+1}={\bar{v}}^{n}+\ldots +\frac{\delta }{\delta y}({\nu }_{d}{\nabla }^{2N}D)+\frac{\delta }{\delta y}({\nu }_{2d}D)\end{array}\right.,$$
(23)

where \(D\) is the horizontal divergence along a Lagrangian surface. The divergence damping consists of the second-order damping \(\frac{\delta {(\nu }_{2d}D)}{\delta x,\delta y}\) and high-order damping \(\scriptstyle\frac{\delta \left({\nu }_{d}{\nabla }^{2N}D\right)}{\delta x, \,\delta y}\) with \(N=2\) used in the operational HAFS. The damping coefficient of the second-order damping is parameterized as \({\nu }_{2d}={c}_{s}A\Delta t\sqrt{{D}^{2}+{\varsigma }^{2}}\), which is similar to the 2D Smagorinsky model35 but with the unit of \({m}^{2}\). \({c}_{s}\) is the Smagorinsky constant, \(\varsigma\) is the vertical relative vorticity, and \(\Delta t\) is the time step of model integration. The damping coefficient of the high-order damping is parameterized as \({\nu }_{d}={({d}_{2N}{\rm{A}})}^{N+1}\)20 with the unit of \(\scriptstyle{({m}^{2})}^{N+1}\), where \({d}_{2N}=0.15\) is an empirical coefficient. It is evident that \({\nu }_{d}{\nabla }^{2N}D\) and \({\nu }_{2d}D\) have the same unit of \({m}^{2}{s}^{-1}\). Apparently, the horizontal turbulence diffusion parameterized by the divergence damping in Eq. (23) is independent from the vertical turbulence diffusion parameterization based on the predicted TKE.

From the perspective of turbulence parameterization, both divergence damping and local down-gradient K-closure are the means to parameterize the unclosed turbulent fluxes induced by the unresolved eddies. Thus, a comparison between them may shed new light on developing consistent horizontal and vertical turbulent mixing parameterization within the framework of FV3. Using the local down-gradient transport formula, for \(i,{j}=1,\,2\), Eq. (1) may be rewritten in the same format as Eq. (23) to yield (Supplementary Note 7),

$$\left\{\begin{array}{ll}{\bar{u}}^{n+1}={\bar{u}}^{n}+\cdots +\Delta t\left[{k}_{m}^{(h)}{\nabla }^{2}\bar{u}+\frac{1}{2}\frac{\delta {k}_{m}^{(h)}}{\delta x}\left({S}_{11}-{S}_{22}\right)+\frac{\delta {k}_{m}^{(h)}}{\delta y}{S}_{12}\right]+\frac{\delta \left(\Delta t{k}_{m}^{(h)}D\right)}{\delta x}\\ {\bar{v}}^{n+1}={\bar{v}}^{n}+\cdots +\Delta t\left[{{k}_{m}^{(h)}\nabla }^{2}\bar{v}+\frac{1}{2}\frac{\delta {k}_{m}^{(h)}}{\delta y}\left({S}_{22}-{S}_{11}\right)+\frac{\delta {k}_{m}^{(h)}}{\delta x}{S}_{21}\right]+\frac{\delta \left(\Delta t{k}_{m}^{(h)}D\right)}{\delta y}\end{array}\right.,$$
(24)

The last term in Eq. (24) is equivalent to the second-order divergence damping with \(\Delta t{k}_{m}^{(h)}\propto {\nu }_{2d}\). This provides a way to parameterize the damping coefficient \({\nu }_{2d}\) in terms of the predicted TKE consistent with the vertical turbulent mixing parameterization. Given the grid spacing of operational HAFS, for the HAFS experiments with the 3D-SA TKE scheme performed in this study, we parameterize \({\nu }_{2d}\) in Eq. (23) in terms of TKE via \({k}_{m}^{(h)}\). Meanwhile, we leave the high-order damping terms untouched. In the future, we will explore venues for applying TKE-based horizontal diffusion to the full momentum equation and other physical quantities (potential temperature, air mass, etc.).

Finally, the implementation of the 3D-SA TKE scheme in HAFS involves both FV3 and CCPP. Specifically, the components involving 3D calculations are implemented in FV3, whereas blending the 3D scheme with 1D scheme takes place in CCPP. Global variables are defined to link FV3 with CCPP. A flow chart of how the 3D-SA TKE is implemented in FV3 and CCPP is shown in Supplementary Fig. 6.