Reflections on the Scientific Legacy of Sergej S. Zilitinkevich on PBL and Urban Parameterizations in NWP Models

The paper summarizes many of the scientific achievements of Professor Sergej S. Zilitinkevich (1936–2021). It first focuses on his basic and applied atmospheric boundary layer research contributions. It then reviews their applications within research and operational numerical weather prediction and air quality modeling, showing their contribution to solving modeling problems related to extremely-stable and -unstable boundary layers.

In 1990 Sergej moved to Western Europe, first as visiting professor at RISOE Nat-ional Laboratory in Denmark. From 1991 to1997 he then had visiting professorships at the Max Planck Institute for Meteorology at the University of Hamburg; Alfred Wegener Institute for Polar and Marine Research; and the Institute for Hydrophysics at GKSS Research Centre in Geessthacht. From 1998 to 2003 he was professor and chair of meteorology at Uppsala University, and since 2004 was a research professor at the Finnish Meteorological Institute, Marie Curie Chair and later research director of European Union projects at the Institute of Atmospheric and Earth System Research at the University of Helsinki. Since 2012, Sergej was also chief scientist of the Pan-Eurasian Experiment, which involves groups from Western Europe, Russia, China, and Japan with a focus on climate change and global pollution. Sergej was an invited plenary speaker at the 10th International Conference on Urban Climate (ICUC10) in 2018. His research contributions have been honored by many international prizes and awards, including the: 2000 Wilhelm Bjerknes award from the European Geophysical Society, 2015 Alfred Wegener award from the European Geophysical Union, and the 2019 International Meteorological Organization (IMO) Prize from the World Meteorological Organization.
The research contributions of Sergej built upon the achievements from the Russian, German, and Nordic schools of geophysical boundary layers and turbulence, e.g., from A. N. Kolmogorov, O. Ladyzhenskaya, A. S. Monin, A. M. Obukhov, L. Prandtl, T. von Karman, V. W. Ekman, C. G. Rossby. His work focused on atmospheric stratified turbulence in both convective and stable layers, with contributions in areas around temperature fluctuations (θ'θ') and buoyancy fluxes. These were latter developed by him into theoretical derivations of the: MOST stability functions, potential turbulence energy concept, and energy-flux balance turbulence closures. This theoretical work contributed to the contemporary nature of these fields and have helped improve PBL parametrizations in numerical weather-prediction (NWP), climate, and air quality (AQ) models worldwide.
In 2007 the journal Boundary Layer Meteorology published a special issue entitled "Atmospheric Boundary Layers: Nature, Theory, and Application to Environmental Mod-elling and Security" dedicated to his 70th Birthday. The introductory paper to the volume by Djolov (2007) summarized scientific achievements of Sergej. It concluded that his new theoretical ideas gave rise to numerous developments in NWP, climate, and AQ modeling in urban environments, including optimal energy usage in variable climate and weather condit-ions, as well as renewable energy planning and management. The volume was extended and published as a book in English and Russian (Baklanov and Grisogono 2007) and it was summarized in BAMS (Baklanov et al. 2011). Many collaborators and former students in the virtual research network of Sergej had papers in that book. The final book by Sergej (Zilitinkevich 2013)  The current paper thus provides the perspectives of its authors as to how the scientific achievements of Sergej in PBL theory are now ever more frequently used in research and operational numerical models for weather prediction, AQ, and climate change predictions for urban and rural communities on a variety of scales. It thus first analyzes his scientific legacy in terms of his contributions, and then provides an analysis of the applications of this work on past, current, and future PBL research and modeling. The possible future applications of his work are partly based on personal interactions between Sergej and the authors of this paper.
As the paper aims to reach both students and seasoned researchers, it first presents the basic state of the art in terms of theoretical concepts before the work of Sergej to better appreciate how his work solved problems and thus improved our understanding of turbulence and PBL structure. It then does the same for research and operational NWP models, most of which are based on the Reynolds Averaging Navier Stokes (RANS) equations. The research areas reviewed for this second section mostly address the origins of problems associated with modeling extreme atmospheric static-stability and -instability, as well as those associated with urban atmospheric parameterizations. Some of these problems still exist in current mesoscale models. An excellent introduction to such models is given in the text of Pielke (2013). One of the most widely used RANS models is the Advanced Research version of the Weather Research and Forecasting model (WRF-ARW) described by Skamarock et al. (2019), and some results from this model are used to illustrate ideas discussed herein. Given the vastness of both fields, the discussions depend heavily on the research experiences of the two authors. The paper uses his own pet nickname "SSZ" when referring to Sergej, but not when citing specific journal papers.

Contribution of Sergej Zilitinkevich to Turbulence Theory
Among the many important scientific achievements in atmospheric turbulence theory of SSZ, this section reviews a selection of his contributions to: (i) analytical PBL height formulations, (ii) newly recognized long-lived stable PBL types; (iii) new understanding of 2-D turbulence at high Richardson numbers (Ri), (iv) a new paradigm for non-local transport via coherent convective flow structures: order out of chaos, (v) the New Energy and Flux Budget (EFB) theory, and (vi) applications to urban PBL structure.

Analytical PBL Height Formulations
The upper boundary of the PBL acts as a barrier to prevent (or strongly reducing) vertical dispersion in most of AQ models, and it also influences extreme weather events and microclimates in NWP and climate models, respectively. Calculation of the PBL height (h P BL ) has been the subject of numerous theoretical and experimental studies since Ekman (1905), who identified the PBL as the steady-state boundary layer in a rotating fluid. He derived the rela-tion h P BL ∼ √ K ma / f , where K ma is a reference eddy viscosity and f the Coriolis parameter. The Rossby and Montgomery (1935) formulation h P BL ∼ u * / f , where u * is the friction scale-velocity, is still used as a rough approximation for PBL height (and its other features). This parameter, however, also strongly depends on static stability, traditionally characterized by the vertical turbulent flux of buoyancy (F B ) at the Earth surface, given by: where g is the acceleration due to gravity, T is temperature, and F h and F qv are the vertical turbulent fluxes of heat and specific humidity (q v ), respectively. Three basic PBL types are thus distinguished as: stable, with F B < 0 reducing turbulence; neutral, with F B = 0; and unstable (or convective), with F B > 0 enhancing turbulence. Zubov (1945) showed that the height of the convective PBL (h C P BL ) will grow if F B remains positive, i.e., More general formulations of CPBL height that account for turbulent entrainment at its upper boundary can be found in Zilitinkevich (1991Zilitinkevich ( , 2012. For stable nocturnal PBLs (SPBLs) over land due to F B changing from positive to negative, Zilitinkevich (1972) derived: h S P B L = C S P B L u 2 * (F B f ) −1/2 , where C S P B L = 0.5 is a dimensionless empirical constant.

Newly Recognized Long-Lived Stable PBL Types
Until recently, atmospheric PBLs were classified as stable, neutral, or unstable, in accordance to near-surface stratification without regard to processes its upper boundary. This is not surprising, as boundary-layer meteorology basically focuses on mid-latitudinal PBLs over land, characterized by large diurnal variations, with typically unstable stratification during daytime and stable stratification at night. Typical nocturnal stable PBLs develop against the well-mixed, and thus almost neutrally stratified residual layer that remains from the daytime convective PBL (CPBL). This scheme is presented in all textbooks, with the result that "nocturnal" and "stable" SPBL are often considered synonyms.
These simple categories were, however, shown as more complex by Zilitinkevich et al. (2006Zilitinkevich et al. ( , 2013, who realized that diurnal variations of surface-T are weak during long polar nights and over the open sea. These periods can be several days or weeks over oceans and can extend to more than a month over both oceans and land during the polar night or day. During such periods, residual layers do not occur, and the PBL then develops against the free tropo-sphere, where global horizontal heat transfer permanently re-establishes stable stratification with a typical Brunt-Vaisala frequency (N) of~10 -2 s −1 (Stone 1972(Stone , 1974 Interactions between the PBL and free troposphere thus suppresses turbulence in the upper PBL and dramatically reduces its height . Long-lived PBLs with no buoyancy flux at the surface were traditionally considered as neutral, but they are in fact strongly affected by the background stratification in the troposphere, which makes them several times shallower than truly neutral PBLs (Zilitinkevich and Esau 2002). Similar long-lived stable PBLs are also several times shallower than usual nocturnal stable PBLs (Zilitinkevich et al. 2007a). This shallowing (Fig. 1) shows atmospheric and LES PHB height values across three stable PBL types: normal nocturnal, marine, and polar. Results show systematic PBL shallowing with increasing values of stability, as parameterized by N. Zilitinkevich and Esau (2002) and Zilitinkevich et al. (2007a) showed that the major forcings and properties (e.g., depth) differ for the two PBL types, the: (i) short-lived mid-and low-latitude continental PBLs have pronounced diurnal variations and thus are almost fully controlled by F B and (ii) long-lived high-latitudinal PBLs and those that keep their stratifycation (stable or unstable) are essentially controlled by the persistent stable stratification in the free troposphere. They are characterized by N~10 -2 s −1 ), making them shallow and hence sensitive to external impacts, such as larger scale phenomena.
Parameterizations of stable PBL processes is a considerable source of uncertainty in NWP, climate, and AQ modeling. The latest models have high vertical resolutions within the PBL and require good surface boundary layer (SBL) parameterizations, where turbulent fluxes are assumed invariant with height. Such parameterizations are especially needed in urban-scale AQ models, where low-level emissions (e.g., from traffic) within the SBL are typical. Recent increases in computer power allows for increasing model grid resolutions, producing better representation of PBL processes, which requires an increased knowledge of the physical, chemical, and biological mechanisms controlling these processes. Normalized stable PBL height (h e ) versus normalized Brunt-Väisälä frequency (N) across three types of stable PBLs. The data are compilation of observations (blue) and LES (red) results, and the horizontal blue line represents the division of applicability between the traditional theory and that of Zilitinkevich et al. (2007a) shown by the solid black line

PBL Turbulence at High Values of Ri
The free troposphere above the PBL is characterized by weak velocity shears and stable stratification, which strongly reduce turbulent intensities. This strongly reduces vertical exchange processes at the PBL upper boundary. Impacts from the Earth surface (e.g., warming or cooling, stronger or weaker surface drag, pollutant emissions) thus almost immediately affect the entire PBL, but only slowly penetrate the free troposphere, especially with unstably stratified, strongly mixed PBLs. The resulting boundary can be seen clearly in vertical profiles of mean T and pollutants from the surface, which almost completely remain within the PBL.
Many diffusion closures include only down-gradient transport, implying that turbulent fluxes are proportional to, and oriented along, the mean gradients of transported properties. The proportionality between the eddy mixing coefficients, i.e., viscosity (K M ), conductivity (K H ), and (pollutant) diffusion (K D ), are the major unknowns in turbulence closure theory. Current closures are usually based on paradigms of Kolmogorov (1941aKolmogorov ( , b, 1942 for the: (i) budget equation for turbulent kinetic energy (TKE) per unit mass (E K ) and (ii) definition of a turbulent velocity scale (u t ,) as equal to E K ½ . The Prandtl vision for neutrally stratified boundary-layer flows was of equal turbulent exchange coefficients (K M~K H~K D ) proportional to the productu t l t , where l t is the turbulent-length scale. Parameter l t was assumed as a linear function of height (z) above the surface, while u t was considered as invariant with z. This relationship has resolved the problem of TKE dissipation rate (ε K ) to be given as: The formulation thus yielded a closure for neutrally stratified sheared flows that revolutionized the understanding and modeling of their turbulence. The followers of Kolmogorov extended the theory (without proof) to stably and unstably stratified flows, giving rise to a turbulence closure theory that is still used in PBL modeling, and in moderately stable or unstable stratifications the approach yields good results. It fails, however, in strongly stable stratification, in which it erroneously prescribes a degeneration of turbulence in the supercritical stratifications typical of free troposphere and ocean thermoclines. Conventional theory and closure models also become erroneous when applied to horizontal diffusion in unstable stratification, such as with pollutant plume spread. A demonstration of the failure (Fig. 2) of conventional theory is its prediction that ε K in stable stratification increases with increasing Monin-Obukhov length (L) via z/L, while the new theory of Zilitinkevich et al. (2007bZilitinkevich et al. ( , 2013 shows that it decreases. In non-neutral flows, TKE production or destruction depends on the gradient Richardson number (Ri), i.e., a non-dimensional stability parameter given as the ratio of the vertical gradient of buoyancy divided by the squared gradient of the mean horizontal velocity (V) sheer: where g is the acceleration due to gravity, θ a a reference potential temperature, and where Ri is positive for stable flows, zero for neutral flows, and negative for unstable flows. For stably stratified flows, the conventional theory predicts turbulence production as possible only in sub-critical stratifications, i.e., when Ri does not exceed its critical value (Ri c ) of 0.25, as obtained by laboratory experimentation. In the atmosphere and oceans, however, sub-critical stratification is observed only with strong velocity shears in boundary layers and jet flows. For supercritical stratifications, the conventional theory predicts a degeneration of turbulence and thus laminarization of the flow. This conclusion, however, contradicts experimental evidence that turbulence is present in the free atmosphere and deep ocean, even with extremely strong supercritical stratifications, i.e., up to Ri~20. For the small-scale flows typical of lab experiments, Ri c > 0.25 is, however, the correct criterion for laminarization. Over   almost a century, this paradox remained unsolved; and the supercritical stratified tur-bulence inherent in the atmosphere and oceans was observed without knowledge of its nature. Zilitinkevich et al. (2013) explained this paradox by illuminating the differences between the small-scale (i.e., low Reynolds number, Re) flows in laboratory experiments and the largescale (high Re) flows in atmospheric/oceanic circulations. Their idea was that turbulence is extinguished only by viscosity, so that the correct criterion of laminarization is based on the turbulent Reynolds number (Re t ). This explains why a supercritical stratification, Ri > Ri c , is sufficient to reduce the already small Re t below its critical value and thus extinguish turbu-lence in low-Re flows, but that it is insufficient in high Re flows. Resulting nocturnal turbulence (q 2 ) is shown in Fig. 3 from Zilitinkevich et al. (2019) for the: (i) traditional 3-D q 2 that decreases as Ri approaches 0.2 and (ii) new weak 2-D horizonal wave-like q 2 that exists with Ri up to about 10.

Paradigm for Non-local Transport via Coherent Convective Flow Structures: Order Out of Chaos
Conventional turbulence theories were based on the paradigm from Kolmogorov (1941a, b), which comprises the following universally recognized postulates: (i) Turbulent flows consist of two types of motion: (i) fully organized mean flows and fully turbulent eddies generated by mean flow instabilities, (ii) eddies are unstable and break down into smaller eddies, thus cau-sing a forward cascade of TKE and other properties from larger to smaller scales towards their viscous dissipation, and (iii) turbulence transfers of momentum, temperature, and air constituents is along their mean gradient, so that their turbulent fluxes are defined as the product of their gradient times a turbulent exchange coefficient, similarly to molecular transport. Zilitinkevich (1973) revealed that theories based on this paradigm, e.g., Monin-Obukhov Similarity Theory (MOST), yield a false picture of turbulence in unstable stratification, but at the time the paradigm was not questioned, and thus the finding was ignored. Zilitinkevich (2010) and Zilitinkevich et al. (2013Zilitinkevich et al. ( , 2019Zilitinkevich et al. ( , 2021, however, theoretically demonstrated why (and could show experimentally that) contrary to postulates (i) and (ii), buoyancy-generated chaotic vertical plumes fundamentally differ from shear generated eddies. The latter are indeed dynamically unstable and break down into smaller eddies, thus performing a forward cascade of 3-D TKE towards its viscous dissipation. Plumes, however, merge to yield larger plumes, thus performing an inverse cascade from smaller to larger scales.
The inverse cascade of buoyancy generated TKE converts it into the mean kinetic energy of large-scale self-organized structures, thus feeding the well-known cells and rolls typical of convective boundary layers, as well as of tropospheric jet streams. A summary of this new paradigm of unstably stratified turbulence (Zlitinkevich et al. 2021) was brought to publication posthumously by four of his Ph.D. students/co-authors. It elucidates the inconsistencies of the conventional theory and explains why that theory still could, however, correctly predict buoyancy generated TKE structure and its vertical exchange coefficients, i.e., because the rate of viscous dissipation of buoyancy generated TKE of convective origin from the traditional interpretation is mathematically indistinguishable from its rate of conversion into the energy of self-organized structures in the new understanding. Zilitinkevich et al. (2006Zilitinkevich et al. ( , 2007bZilitinkevich et al. ( , 2008aZilitinkevich et al. ( , 2009Zilitinkevich et al. ( , 2013Zilitinkevich et al. ( , 2019 created the non-conventional Energy and Flux Budget (EFB) theory of stably stratified, high-Re (geophysical/astrophysical) turbulence that explains both subcritical and supercritical stratifications. Its key ideas from the counter-gradient mechanism of temperature and buoyancy transfer [disregarded in Postulate-(iii)] and from feedbacks from this mechanism, show that: (i) a negative buoyancy flux does not extinguish TKE, as commonly believed, but converts it into turbulent potential energy (TPE), proportional to the mean-squared temperature fluctuation; (ii) warm eddies (positive temperature fluctuations) acquire positive buoyant acceleration and rise, while cold eddies (negative fluctuations) acquire negative acceleration and sink, so that both transfer temper-ature upward, against its basic downgradient flux, and (iii) this yields a self-limitation of the temperature and buoyancy fluxes, assuring the maintenance of turbulence in supercritical stratification.

New Energy and Flux Budget (EFB) Theory
The EFB theory illustrates the difference between the Prandtl number (Pr), given by K m /K H , in the two types of stably stratified turbulence, i.e., in the: (i) familiar subcritical regime, the overall-mixing turbulence has similar vertical exchange coefficients for momentum and heat, so that Pr is nearly constant and (ii) in the supercritical regime, the newly discover-ed weak-conductivity turbulence has heat transfer and diffusion strongly damped by stratification, so that Pr~Ri. Atmospheric and laboratory measurements, as well large eddy simulation (LES) and direct numerical simulation (DNS), estimates of Pr values in both the traditional and wave-like stable-stability regimes are shown in Fig. 4 (from Zilitinkevich et al. 2013). Results show the expected value of unity in the traditional regime and increasing values in the new regime. Meteorologists and oceanographers have long known that turbulent heat and mass transfer is weaker than momentum transfer in the free atmosphere and deep ocean. The EFB theory, based on the insight that the total turbulent energy (sum of the TKE and TPE, proportional to the potential temperature variance) thus explains, quantifies, and offers a physically grounded method for modelling high-Re geophysical turbulence in any stable stratification.
In summary, the work of SSZ has made it possible that atmospheric turbulence can now be resolved for stable stratification via the EFB turbulence closure. The problem of unstable  Zilitinkevich et al. (2013) stratification can now likewise be resolved via his concept of self-organization of convective turbulence.

Urban PBL Structures
Urban weather, AQ, and climate models have a variety of parameterizations associated with urbanization effects on different scales. They also have operational requirements, such as forecasting versus assessment, or environmental AQ applications versus emergency preparedness. Since the end of 1990s SSZ was actively involved in studies of urban boundary layers (UBLs), in particular with European projects like FUMAPEX (Baklanov et al. 2008) and MEGAPOLI (Baklanov et al. 2010). Incorporation of urban effects into urban-and regionalscale NWP and AQ models use the urban PBL parametrizations of SSZ, such as urban his: (i) aerodynamic roughness-length (z 0 ) ratios for momentum, heat, and tracers as a function of roughness Reynolds number (Zilitinkevich 1970;FUMAPEX 2003), (ii) effects of stratifica-tion on z 0 and zero-plane displacement height (d 0 , Zilitinkevich et al 2008b), (iii) h SSBL values from the  parameterization in combination with a prognostic advec-tion and diffusion equation , (iv) twoway interactions between UBL turbulence and aerosols Kulmala et al. 2023), (v) imple-mentations in a mesoscale: NWP models, such as DERMA, HIRLAM, Enviro-HIRLAM (Baklanov et al. 2017), WRF; projects, such as FUMAPEX, MEGAPOLI, Pan-Eurasian Experiment (PEEX, Kulmala et al. 2015); World Meteorological Organization (WMO) training and education efforts (WMO 2020), and (vi) WMO studies of personal environments  for new services via the Integrated Urban Systems and Services Program (WMO IUS 2019, 2021). In his last years, SSZ was quite interested in contributing to such efforts. He thus helped the University of Helsinki team to win the Russia MEGAGRANT on urban air quality and made significant contributions to science-policy interfaces for the final TRAKT-project publication (Esau et al. 2021).
To expand the urban related points, Zilitinkevich et al (2008b) demonstrated that urban z 0 and d 0 values characterize the resistance exerted by its surface roughness elements on turbulent flows for a wide range of turbulent-flow problems. Classical laboratory experiments and theories treated both parameters as geometric parameters, and thus as independent of flow characteristics. They demonstrated their stability dependences, stronger for z 0 (especially in stable stratification) and weaker (but still pronounced) for d 0 . They also developed a scaleanalysis model for these dependences and verified it against experimental data for urban and forest canopies. They proposed that stability effects on z 0 can be incorporated via z/L: where ϕ m and ϕ h are non-dimensional SBL stability functions for momentum and heat, respectively, and z/L is a non-dimensional stability parameter that defines the height above the ground at which buoyancy produced q 2 first becomes greater than shear produced q 2 . Both functions are unity in a neutral, equilibrium SBL (i.e., one with equal sources and sinks of TKE). The non-neutral stability equation for z 0 for momentum (z om ), was thus related to average urban building height (h b ) by: where C ss and C us are constants equal to 8.13 and 1.24, respectively. Urbanization amplifies local environmental stresses, which strongly increases damages from natural disasters. National Weather Services (NWSs) provide meteorological information averaged (at best) over 3 km squares. Urban morphologies, however, produce strong small-scale heterogeneities of T, V, and precipitation, and thus peak damages during natural disasters result from local extremes produced by these heterogeneities. The major challenge for NWSs is thus to move from centralized, limited-resolution weather, climate, and AQ services that only address large bulk segments of society towards client-oriented services capable of resolving the fine features of extreme and dangerous events at human and community scales. Zilitinkevich et al. (2015Zilitinkevich et al. ( , 2016 in collaborations with the WMO thus suggested a focus on the physics of personal environments, i.e., the climate, weather, and AQ within 100-m sized atmospheric cubes. It stressed interactions and feedbacks between physical and chemical processes within an UBL, e.g., local pollutant emissions alter its radiation budget, so that it becomes more stable. This suppresses turbulence, hence decreasing it height, and thus yielding even stronger concentrations and even lower PBL heights (Baklanov et al. 2008;Petaja et al. 2016;Ding et al. 2016).
Complex unique climates form over the highly heterogeneous urban surfaces. Their most studied feature, i.e., the urban heat island (UHI), is defined as the T-excess of the city over its surrounding rural areas. Oke et al (2017) lists several UHI types: subsurface, obtained from buried sensors; surface (radiative), obtained from overhead observations of emitted long wave radiation; 2-m air, obtained from a network of T-sensors; and PBL, obtained from soundings. To evaluate mesoscale NWP model output on its horizontal grid spacing ( x) of 1 km with 2-m urban T-data is difficult, as most networks are too sparse, so surface radiative UHIs can be useful if their resolutions are neither too coarse nor detailed.
Validation of NWP output with satellite derived surface temperatures in urban areas, however, is not easy. Under the best circumstances, when the sensor and sun (during daytime hours) are both directly overhead to avoid shadow effects (Voogt and Oke 2003), the surface T-values from the satellite are essentially roof and road values, so it is important to ensure that the NWP output represents the same. Such models thus require a surface scheme able to differentiate between roof, road, and wall values (Stewart et al. 2021). Figure 5 shows a 1992 LANDSAT-5 IR image over Washington D.C. (from S. Stetson, personal communication) at a resolution of 150 m overlayed on an aerial photography at a resolution of 0.33 m. The features of its surface radiative UHI are detailed enough to show the individual buildings forming the personal environments envisioned by SSZ, yet large enough to be validated with output from an NWP model.
Understanding and modelling local urban atmospheric environments requires a multidisciplinary approach that integrates the revised turbulence and PBL theories of SSZ into multidisciplinary research on the physics and biogeochemistry of atmosphere-ecosystemmegacity interactions. The WMO thus foresees that the current centralized top-down monitoring by NWSs will be evermore supplemented with bottom-up observations by individuals, businesses (such as transport, agriculture, energy) and city administrations (Mahoney and O'Sullivan 2013;Grimmond et al. 2014Grimmond et al. , 2015Snik et al. 2014;Baklanov et al. 2018Baklanov et al. , 2020. To support the future use of such crowd-sourced data, a deeper understanding of the processes governing personal environments, as well as new ideas about their modelling and forecasting, are needed (Muller et al. 2015).
An approach that could make use of crowd-sourced data is the application of computational fluid dynamics (CFD) models to urban areas. Such models are run with 3-D grids of 1-2 m, so they can only be applied to neighborhood scales. They have been used to simulate the effects of regional climate changes (Conry et al. 2015) and to test UHI mitigation strategies via white roofs and/or greening (McRae et al. 2020). While neither study assimilated crowdsourced data, nor used them in a model evaluation study, such models are more suited than mesoscale RANS models to do so, given the large amount of data that could be soon available in a single mesoscale model grid cell from mobile phones (Erez Weinroth, MobilePhysics, Inc., personnel communication).

Applications of the Work of Sergej Zilitinkevich in NWP Models
The above section summarized high points of the work of SSZ, and how they have significantly increased our understanding of PBL turbulence characteristics and life cycles. The next section highlights past, current, and possible future applications of this work within the SBL and PBL formulations in NWP and AQ models. It first presents a history/ summary of problems with 1st order turbulence closures in NWP models, first over homogeneous terrain and then over heterogeneous terrain. How the contributions of SSZ helped to address these problems is then reviewed.

First Order Closures
This section first discusses problems associated with specified-shape 1-D vertical profiles for turbulent exchange coefficients, K(z), in early atmospheric models during stable conditions to show the origins of the problems associated with the modeling of strong atmospheric static-stability and -instability. The earliest problems and solutions were reviewed by Shir and Bornstein (1977, herein SB77) and this section borrows from that reference, as some issues are still relevant in models with 1 st order closures. A first PBL model to incorporate 'local stability' effects in K(z) was Stevens (1959), who used model ∂θ/∂z values at z to modify his linearly decreasing (from its SBL value) K(z). As this stability effect produced an unrealistic mid-PBL minimum (discussed below), Estoque (1961) incorporated nonlocal stability effects on his also fixed (but non-linear) shape PBL K-profile via MOST, so that his K-value at his prescribed SBL top (h S ) was: where l m (= kz) is the Prandtl SBL mixing length for momentum and k (= 0.4) is the von Karman constant. The Monin-Obukhov (1954a) stability function used was: where α is a positive "constant" whose value(s) will be discussed below. As SBL stability changed, so did K m (h s ), which altered the PBL K-profile. As the profile excluded local stability influences, the K-minimum mentioned above was avoided. Fixed shape profiles were, however, considered less desirable than those that incorporated local stability effects, so Pandolfo et al. (1963) used (3), but with the Blackadar (1962) PBL profile for l(z). They also used ϕ m (z) formulation for forced (z/L > −1) and unstable free (z/L < −1) convection, via the SBL formulations of Priestley (1959). The problem was not, however, eliminated, so Sasamori (1970) proposed the fixed parabolic K(z) profile of O'Brien (1970), as observed by Elliot (1964). The source of the "problem" was identified by SB77. Their simulated V-profile (Fig. 6) produced a low-level jet (LLJ) at a height of several 100 m after 6-h of cooling. The corresponding K(z) profile (via Pandolfo et al. 1963) shows the local minimum, while their θ-profile shows a lack of cooling above the minimum. The impediment of nocturnal downward heat flows that result from the K-minima produces both a too low surface-T (Sasamori 1970) and urban circulations confined within the layer below the minimum (Bornstein 1975), both still problems in some NWP models.
They further showed that the minimum arises when the Ri is evaluated near a LLJ peak. The associated weak shear produces a too large "finite difference" gradient Ri (i.e., Ri ), given as the following numerical analog to (1): where z (= z 2 −z 1 ) involves V and θ at two adjacent model grid levels. If ϕ m = ϕ h in stable conditions. (2) and (4) Shir and Bornstein (1977) The explanation for the K-minimum is thus: as Ri grows towards 1/α at the LLJ nose, 1-αRi decreases, then (ϕ m) −2 in (3) increases, and finally K decreases. Peak (and not minimum) K-values derived from turbulent flux measurements were, however, observed near LLJ peaks within the elevated west-coast subsidence inversion by Lile (1970).
The maximum allowable stable-PBL value of Ri or Ri in (6), i.e., the "critical" (Ri c ) is thus 1/α, whose correct value in the PBL was uncertain (Panofsky and Prasad 1965). Observations from 16 m SBL flux towers yielded values of α from 0.3 to 7 (Neumann and Mahrer 1971), with a value of α = 5 now generally accepted. This compares well with the value of 0.2 for Ri c seen in the laboratory observations summarized by SSZ, discussed in Sect. 2.3.
A small Ri c was a problem in early NWP models, as their lowest model level was generally greater than the 16 m tall flux towers. Equation (5) shows that Ri increases with layer thickness (Reiter and Lester (1967), and thus in many early K-theory models, SBL and/or PBL Ri -values too quickly exceeding Ri c after sunset. The problem was avoided in a variety of ways, including use of a: (i) α < 5, e.g., Estoque (1961) used 1.0 andMcPherson (1970) used 0.03, (ii) reduced Ri value, e.g., Tag (1969) divided his by five, (iii) new "stability parameters," which replaced the squared shear in Ri by a linear shear (Estoque and Bhumralkar 1968), and (iv) arbitrary maximum value for Ri or a minim-um value for K(z). Such fixes are still used in some NWP models with TKE closures, which commonly have lowest model levels of 10-20 m, although they frequently are not mentioned.
Another approach was use of new SBL (ϕ m ) −2 functions in (3) that did not increase as rapidly in stable conditions as those in Monin-Obukhov (1954b). Some of the newer ones (Fig. 7,   improvement, as they produce realistic analytical SBL fluxes that provide the numerical PBL above with better lower boundary conditions (BCs). Lyons et al. (1964) conclude that Ri values are only valid stability-estimates in shallow layers with strong monatomic shears. The unreliability of calculated Ri values in Fig. 8 show that small variations in the radiosonde V-profile produce multiple wide swings in Ri . McNider and Pielke (1981) thus concluded that α (and hence Ri c ) are not constants in NWP models but are related to their z by the following empirical relationship: Ri c = 1/α = 0.115( z) 0.175 .
For unstable daytime conditions, local Ri values also underestimate turbu-lent mixing, as MOST breaks down at about z/L < -1, as the SBL transitions from (shear) forced-to freeconvection. Various unstable forced convection formulations for ϕ m have been proposed. Paulson (1970), for example, used one that depends on a "bulk" Ri (i.e., Ri b ), in which the z of (5) is replaced by z. During unstable forced convection, the PBL is still dominated by horizontal motion via V and u * , as it was during stable conditions. During free convection, however, it is dominated by vertical motions and by w *, given by Deardorff as w * = (cg F h h C P BL ) 1/3 , where c = 1.2.
During free convection, the PBL is thus well mixed and has a near neutral stability, even with strong local mixing from rising thermals from the super-adia-batic SBL. Such non-local PBL mixing cannot be captured by K(z) formulations dependent on local gradients. The thermals are too large to be modeled by the diffusion term and too small to be captured by the vertical advection term (Ackerman and Appleman 1975). PBL K(z) formulations dependent on SBL gradients thus again better represent eddy processes in convective PBLs. Kuo (1968) showed such K(z) profiles as depended, not on local ∂θ/∂z values, but on the integra-ted stability from the surface to the level where θ equals its surface value. The "penetrative convection" term of Estoque and Bhumralkar (1969) was an early attempt to parameterize non-local convection processes in a 3-D PBL model.
In summary, while standard 1-D K-theory over homogeneous terrain were adequate under neutral conditions, no acceptable formulations existed for strongly non-neutral conditions until the work of SSZ discussed in Sect. 2.4. Additional problems over non-homogeneous terrain are discussed next. . Also shown is the shape of the resulting equilibrium (EBL) and internal (IBL) boundary layers. See text for explanation of regions a-d; from Shir and Bornstein (1977) Results from the 2-D TKE model of SB77 showed that K(z) profiles in a neutral-stability flow in the equilibrium (i.e., equal sources and sinks of TKE) regions (a) and (d) of Fig. 9 depend only on z' 0 and z 0 , respectively. Those in region (c), the nonequilibrium region between the equilibrium boundary layer (EBL) with its slope (i.e., vertical depth divided by distance downwind from the discontinuity) of 1/100 and the internal boundary layer (IBL) with its slope of 1/10, depend on both values. In region (b), above the IBL, the K(z)is related only to the upwind z' 0 and not to z 0 , as effects from the new surface has not yet reached that level. Their results also showed that ϕ m values from (4) are unity only in the equilibrium regions of (a) and (d). Peterson (1971) showed that K(z) profiles over inhomogeneous terrain thus must also depend on horizontal gradients. The EBL slope also implies that for MOST theory to be valid, the ratio x/ z near the surface in NWP models must be < 100 so that h SBL is within the EBL, e.g., for a x of 1000 m, the first model z must be at most 10 m.
In summary, 1st order K-theory is unable to capture vertical PBL structures in strongly stable and strongly unstable conditions over homogeneous terrain. It is also unable to capture them over inhomogen-eous terrain under any stability conditions. This led to application of 2nd order TKE PBL closures in NWP models.

Second Order Closures
A hierarchy of higher order turbulence closure levels for the RANS equations in NWP models was presented by Mellor and Yamada (1974). Not counting the 1 st order prognostic partial differential equations (PDEs) for mean momentum and θ, and ignoring water in all forms, they obtained the 2nd order PDEs and diagnostic algebraic equations (AEs) for different closure levels: (i) Level 4, with no simplifying assumptions, yielded 10 PDEs: six for the independent stress components (τ), three for the heat flux components (F h ), and one for potential temperature variance (θ 2 ); (ii) Level 3, ignoring some "higher order" terms, but keeping advection and horizontal diffusion, yielded two PDEs for θ 2 and TKE (written herein as e) and eight AEs: for the five remaining independent τ-components (after TKE is known), and three for the F h components; (iii) Level 2, also ignoring advection and horizontal diffusion, yielded no PDEs and only 10 AEs; and (iv) Level 1, only 1st order closures. Mellor and Yamada (1982) added a Level 2.5 closure: like Level 3 but replacing the prognostic θ 2 equation with an AE. Some 2.5 closure models also include PDEs for potential turbulent energy (PTE), ε K , l, e/l, and/or h. SSZ thought, however, that l and h were not conservative and thus could not be evaluated by prognostic equations. SSZ further believed that prognostic equations for ε K had to be "constructed" analogous to the TKE-equation, as the physics of each of its source and sink terms was unknown. He thus believed the best NWP models had PDEs for TKE, PTE, and θ 2 , as well as AEs for , l and h.
The 2nd closure e-ε 1-D neutral PBL (NPBL) model of Freedman and Jacobson (2002) found an error in basic e-ε NPBL formulations (e.g., Xu and Taylor 1997) in which the turbulent PBL eddy coefficients were defined as: The constants σ e and σ ε were given values of 1.0 and 1.3, respectively, implying that K e = K m and that K e < K m . Their theoretical analysis, however, showed the values as 1.6 and 1.1, respectively, which alters both relationships to produce: K m > K ε > K e , which thus alters the relative rates of decrease with z of the K-parameters. Their simulations showed that the original l(z) profiles continue to decrease to the model top (Fig. 10a), while the new ones reach zero around that level (Fig. 10b), consistent with the widely used Blackadar (1962) equation and more consistent with the 1-D DNS results of Coleman (1999). The new results imply that mechanical eddies do not continuously grow with z. The new l-profiles produce more reasonable K m -profiles than the old formulation and again better agree with the DNS results (Figs. 10c vs. d). The work was extended to the SPBL by Freedman and Jacobson (2003). Their results showed the same inequality as in the NPBL case, but with slightly altered values for σ e and σ ε . from the earlier NPBL work and those in previous SPBL models (e.g., Apsley and Castro 1997). The two results show that σ e and σ ε are in fact dependent on stability. The work was not extended to the CPBL, due to the non-local nature of its closure.
The National Oceanographic and Atmospheric Administration (NOAA) operational High resolution Rapid response (HRRR) model is based on WRF-ARW. Its MYNN-EDMF PBL option (Olson et al. 2019) uses many ideas from SSZ. The Level 2.5 (and not Level 3) version of the model (hereafter MYNN) was used in a research mode by Melecio-Vázquez (2021, hereafter MV21). He compared its results against those from two other Level 2.5 models with less complex P B formulations, i.e., BouLac (Therry and Lacarrère 1983;Bougeault and Lacarrère 1989) and MYJ (Mellor-Yamada-Janjic) of Janjic (2001). All three models used the WRF MYJ (also called Eta) SBL option (Janjic 2019). That scheme is based on Monin and Obukhov (1954a, b) similarity theory with the Beljaars (1995) correction to avoid singularities in unstable SBLs with vanishing V-values. Their viscous sub-layers over land use variable roughness heights for T and q v , as proposed by Zilitinkevich (1995). MYJ/Eta also imposes an upper limit on the master length scale l (Mesinger 1993). The need to impose an upper limit on l is a necessary replacement of the Mellor-Yamada (1982) limits on the momentum and heat stability functions (Janjić 2019). Otherwise, turbulence generation in cases of intense shear and/or buoyancy production of TKE will be prevented by the division of the TKE generation terms by a too large value of l.
In Level 2.5 (or higher) NWP models. the local rate of change of turbulent energy per unit mass (q 2 /2) results from 3-D advection, 3-D diffusion, shear production, buoyancy production (in unstable conditions) or destruction (in stable conditions) P B , and molecular dissipation (Mellor and Yamada 1974). The focus in MV21 was on closures for P B in the  Coleman (1999). The curves in (a) and (c) are for the old theory, while those in (b) and (d) are for the new theory of Freedman and Jacobson (2002) their three schemes. The P B term in such models can be divided into three processes: where the values of the two constants (β's) differ from model to model and where the subscripts 'up' and 'en' refer to processes (discussed below) within ascending plumes and the environment, respectively. The three processes comprising P B are the: (i) normal downgradient fluxes, as functions of the vertical gradients of θ and q v , (ii) counter-gradient convective buoyancy-flux correction parameterization for heat ( h ) and water vapor ( v ), and (iii) non-local convective mass-flux parameterization (M times final bracket; hereafter M). The last term is included, as it also is a source of q 2 in unstable conditions but is not activated in stable and neutral conditions. In 2nd order closures, PBL eddy diffusivities are no longer functions of mean-quantity gradients, as in (3) for 1st order closures, but are now functions of q, e.g., for TKE: where the TKE mixing length (l e ) and stability function ( Only MYNN has an option for a Level 3 closure and an M-term; while the former was not used in MV21, the latter was. Its e and v formulations were not activated, however, as they are only used when M is not activated. Its M-formulation is based on Soares et al. (2004), which gives it as: M times the difference between the TKE within the ascending plumes ('up') and the environment ('en') in (7). The value of M depends on the fractional horizontal upward plume-area times the w within the plume. Activation requires a positive surface buoyancy flux and a super-adiabatic SNL profile; details are in Olson (2019). The BouLac P B formulation only includes the down-gradient flux of θ (and not of q v ) and thus it only includes H . Its lack of a M-term implies that it can only partially capture the vertical transport of SBL heat in super-adiabatic conditions. This is true, as all K-formulations (TKE or mean grad-ient based) can only capture the effects of down-gradient random turbulent diffusion, while non-local M-formulations simulate the upward transport from the self-organizing vertical convective eddies envisioned by SSZ. The BouLac formulation for H is given by Højstrup (1982) and Troen and Mahrt (1986) as: Counter gradient heat fluxes are important when rising convective plumes impinge on a PBL caping inversion base (i.e., with subsidence or the tropopause). This produces downward mixing of warm air into the PBL in the "entrainment zone" layer, which encompasses the inversion base. Counter gradient momentum fluxes are also important in sustaining generalcirculation jet streams (Lau et al. 1978). Based on their respective P B formulations in MV21, BouLac should be least effective in vertically distribut-ing excess surface heat and thus should have the strongest super-adiabatic layer. MYNN should be most effective, given its M parameterization and it should thus have the weakest super-adiabatic layers.
MV21 model results on a 1 km horizontal grid were tested against data from a day during the Long Island Sound Tropospheric Ozone Study (LISTOS) that produced well defined sea breeze fronts (SBFs) in the New York City (NYC) area. Results from the three PBL models are shown in 1200 LST north-south vertical T-cross sections through "rural" Long Island (Fig. 11), about 100 km east (and not downwind) of the center of NYC. By this time a local SBF had penetrated northward to about the same point near the northern edge of the Island in all models. MYNN shows the strongest upward motions at the front, while BouLac has the weakest, with MYJ as intermediate. As anticipated, MYNN shows the weak-est super-adiabatic SBL, while BouLac has the strongest, with MYJ again intermediate. Error Limits exist for RANS TKE models, e.g., the "terra incognita" regime identified by Wyngaard (2004) who showed that as x decreased in NWP models, the separation between the mean and turbulent flows inherent in Reynolds averaging fails. This occurs at about x = 400 m, where the most energetic (turbulent) eddies are neither fully parameterized in traditional RANS models nor fully resolved in traditional LES models. SSZ believed that RANS Fig. 11 continued with x values near the critical value must use a full 3-D Mellor-Yamada Level 3 closures with AEs for the divergence of the off-diagonal horizontal fluxes (in addition to prognostic equations for TKE and θ 2 ) to obtain reasonable results. This idea was tested in idealized WRF simulations by Juliano et al. (2021). Results were better than when the AE terms were parameterized by K-theory, and he also high-lighted the need for improved l-formulations, an area pioneered by SSZ (Perov et al. 2001;Woetmann 2010).

PBL Height Determination from NWP Output
Diagnostic PBL height is an important product in NWP models for many reasons, e.g., it is the height to which pollutants can vertically spread by diffusion. Various methods are used for this diagnosis, e.g., the height of the nocturnal stable PBL (h SPBL ) is frequently set as the top of the surface radiation inversion, as it the level where surface TKE first decreases to either a prescribed minimum or fixed percentage of its surface value. The height of the convective PBL (h CPBL ) is frequently found from predicted mean T-profiles or (again) from a prescribed critical TKE magnitude.
The 1-D e-l model of Freedman (1996) simulated nocturnal PBL growth, with results tested against data from day 33 of the Wangara (Australia) field data (Dyer and Hicks 1970). Results showed that h S P B L , estimated as the level where TKE was reduced to 5% of its surface level, was not at the surface inversion top but was only about half that height (Fig. 12a). This occurred as turbulent cooling only dominated in the lower half of the inversion, while radiative flux divergence dominated in its upper half (Fig. 12b). Methods that do not use TKE values to determine z SPBL height, as the top of the nocturnal surface inversion can only thus underestimate this parameter.
They also showed the importance of the elevated nocturnal residual layer (RL), i.e., a layer of left over daytime TKE above the h S P B L (Fig. 13). Their TKE profiles show maximum  Freedman (1996) values of about 10 m 2 s −2 during daytime, 0.1 m 2 s −2 near the surface at night, and 0.01 m 2 s −2 in the RL. The peak TKE in the RL reduces with time, while its upper and lower limits converge, eventually only producing a series of sporadic TKE bursts centered at its midpoint elevation. They also showed that the K D -values (not shown) associated with the new TKE and l values resulted in diffusion coefficients comparable to those of the stable Pasquill-Gifford-Turner sigma-curves (Venkatram 1996) used in many AQ models. It is thus necessary to  Freedman (1996) develop methodologies to also determine the heights of the upper and lower levels of the nocturnal RL. Implementation of SSZ suggested parameterizations for SPBL height (see Sect. 2) into AQ models provided substantial improvements of forecasts in extremely stably stratified PBLs . LeMone et al. (2013) evaluated eight objective methods in four PBL schemes available in WRF to determine h C P BL values from criteria based on thresholds of θ, Ri, or TKE. Results were compared with PBL observations from the 1997 Cooperative Atmosphere-Surface Ex-change Study (CASES-97) with its extensive network of rural observations during dry con-vective periods. They tested predictions of three possible estimates of h C P BL height (Fig. 14), i.e., the: bottom (h 1 ) or top (h 2 ,) of the entrainment layer (ET), and minimum (negative) normalized heat flux (z i ,) near the center of ET; TKE profiles do not have a region of negative values necessary to determine z i . The criteria used either a: parcel method (discussed above), maximum Ri value, or minimum TKE value; results showed all methods with only limited success. LeMone et al. (2014) used data from the same site, but during three moderately windy fair-weather nights to evaluate eight methods in the same PBL schemes to determine stable nocturnal h SPBL values. The method was based on four model profiles of TKE and two of Ri, as well as two subjective methods with TKE profile. Results were com-pared to the level where the observed TKE decreased to 5% of its PBL maximum; results again showed all methods with only limited success.
Conclusions concerning the results reported in this section include that: (i) T-criteria to determine z C P BL is more difficult in the absence of a capping subsidence inversion, (ii) a robust method to determine all PBL heights should give the same results if the search was upward (as currently done) or downward directed, (iii) impediments to success for the second Fig. 14 Vertical profiles of: a virtual potential temperature (θv, K) and b normalized vertical heat flux, and how their estimates of daytime PBL heights (i.e., h 1 , h 2 , and z i ). See text for explanation of symbols; from LeMone et al. (2013) recommendation are slight variations in all observed (see Fig. 8) and some model profiles (see figures in LeMone papers), (iv) a way to minimize this last effect is to smooth the observed or simulated profiles, e.g., following the method of Hägeli et al. (2000), (v) given the large z values in the upper levels of NWP models, the application of (ii) will yield large differences in CPBL heights regardless of search direction, (vi) interpolation between adjacent levels can minimize this last effect, and (vii) standards are needed for the critical values of Ri, TKE, etc. used to determine both CPBL and SPBL heights.

Urban PBL 3-D Structures
Good h SPBL determination is especially important in cities, where urban processes can produce one or more weak, shallow elevated inversion layers (Bornstein 1968;Baklanov 2002). Figure 15 illustrates this for an early morning summer day over NYC, where the T and SO 2 data were collected by instrumented helicopter. The local (low resolution) NWS radiosonde on the southern urban edge showed its lowest inversion base at 1000 m, while the helicopter sound-ings showed urban induced inversions at 300 and 400 m. The lower one formed an impene-trable lid on the vertical diffusion of pollutants emitted near the surface, while a power plant plume penetrated only the lower inversion and was thus trapped between it and the upper one. NWP models must have enough urban physics to capture the radiative and turbulent processes producing such weak urban induce inversion layers.
Urban PBLs (UPBLs) are more complex than rural ones Oke et al. (2017), During daytime (Fig. 16a) the warm /turbulent/polluted mesoscale urban boundary layer (UBL) consists of a surface layer (SL) and mixed (sometimes "mixing") layer. Acronyms and symbols from the figure are used here for convenience, even if they do not coincide with those in other parts of the paper. With a mean wind the UBL forms an advected plume over a downwind rural boundary layer, while without a mean wind an urban dome shaped PBL forms instead (Bornstein 1968). On the local scale, the SL is comprised of the urban canopy layer (UCL, with a top at Z H ) between its buildings, which forms the lower half of the roughness sublayer (RSL, with a top at Z r ). The upper most sublayer of the SL is the inertial sublayer (ISL, with a top at Z L ). The UBL is frequently capped by an entrainment zone (EZ), above which is the free atmosphere (FA). The EZ is frequently centered on the base of an elevated capping inversion (CI, not shown in Fig. 16a). At night (Fig. 16b), the upwind rural nocturnal surface Fig. 15 Early morning (0708-0742 LST) north to south (left to right) vertical cross section (generally along mean wind, V) over New York City, with SO 2 concentrations (ppm) obtained from an instrumented helicopter. Cross hatched areas represent elevated inversion layers (I and II), as determined from the helicopter temperature data. Vertical line is location of a slight shift in the plane of the section; from Bornstein (1987) inversion layer (NSL) is advected over the city, producing a weak shallow elevated inversion (first observed by Bornstein 1968) and sequentially capped by the nocturnal RL, sometimes a CI, and then the FA.
While the urban UCL can be simulated via CFD, LES, or DNS simulations for a single canyon up to a full neighborhood, the small grids required (1 m for the first two and I mm for the latter) negate these approaches for NWP models. As urban UCLs and RSLs are not in TKE equilibrium, MOST theories are not valid (Oke et al. 2017). Their analytical parameter-izations in NWP models are thus modified rural SBL formulations developed from observa-tions in urban areas, many based on the work of SSZ (as discussed in Sect. 2.3).
Urban NWP and AQ modeling within the European FUMAPEX and MEGAPOLI projects (Baklanov et al. 2008(Baklanov et al. , 2017 utilized PBL improvements from SSZ (see Sect. 2.6). They concluded that improved forecast accuracy requires parameterizations of: (i) assimilated surface characteristics from satellite data, (ii) high-resolution land-use classifications and algorithms for roughness parameters based on morphologic methods, (iii) model downscaling, including increased vertical and horizontal resolutions and (one-or two-way) nesting, (iv) UCL meteorological fields, and (v) urban turbulent fluxes.
The urbanization options in WRF are summarized in Skamarock et al. (2019), they include the: (i) default Noah land use scheme (Liu et al. 2006), which represents urban areas as only changes in surface parameters, such as z o , specific heat, and albedo, as first done in Bornstein (1975). This provides reasonable HUI estimates, but cannot reproduce building drag effects on PBL winds, (ii) Single Layer Urban Canopy Model (SLUCM, Kusaka et al. 2001;Kusaka and Kimura 2004), which solves complex surface energy balance equations within a single urban canyon, and (iii) multi-level PBL Building Environment Parameterization (BEP) of Martilli et al. 2002).   Oke et al. (2017) As buildings in mega-metropolises like NYC can reach hundreds of meters, NWP models need to address the building interaction at multiple model levels. The BEP "porous building" interactions use of known building fractions in each numerical cell. It can be combined with the Building Energy Model (BEM) option (Salamanca and Martilli 2010), which parameterizes building momentum and energy balances. For each NWP model PDE (i.e., momentum, heat, and TKE), BEP adds three new terms to parameterize impacts from streets, roofs, and walls. The first two terms are functions of u * , while the last is a function of building drag coefficient (C d ). As BEP + BEM can overestimate both 2-m UHIs and wind speeds, Gutierrez et al. (2015a, b) added new variable C d and cooling tower (CT) parameter-izations, the latter of which divides anthropogenic heat (Q A ) from commercial buildings into its latent and sensible heat components. The variable C d reduced the over-estimated V-values, while the CT formulation reduced overestimated NYC UHIs.
The urban morphological data for each of the three urbanization schemes (SLUCM, BEP, and BEP + BEM) can be applied to the one, three, or 10 urban classes in the Noah, BEP, and local climate zone (LCZ, Stewart and Oke 2012) schemes, respectively. LCZ data can be obtained from: (i) WRF lookup tables, (ii) the World Urban Database and Access Portal Tools (WUDAPT) methodology Bechtel et al. 2019) that uses remotely sensed surface data to produce grid specific distributions of some required parameters, or (iii) building specific data for use in WUDAPT, e.g., the NYC Primary Land Use Tax Output (PLUTO) data in Ortiz et al. (2018). Global LCZs maps by Demuzere (2022) can support earth system modeling and urban scale environmental science.
While use of BEP + BEM can double WRF CPU times, the results of Martilli et al. (2002) show its value. They used the BEP and Noah urbanization schemes and compared both with urban canyon TKE observations (Fig. 17). The data were a synthesis by Rotach (2001) of observations over Zurich, Basel, and Sapporo, and of wind tunnel results. The Noah results showed expected TKE decreases with height during both day and night hours, with the latter values reduced over the former. The new results showed TKE: (i) reaching a maximum just above the grid level with the maximum number of buildings, (ii) reduced within the UCL and increased above rooftops, as compared to Noah values, (iii) with identical profile day and night, showing dominance of mechanical over thermal TKE production, and (iv) with a good match with the observations.
The BEP + BEM + variable C d + CT urbanization package (BEP + BEM) + , hereafter], not yet in NCAR WRF releases, was used in the NYC LISTOS summer day study of MV21 (discus-sed above). The package generally produced more accurate urban T and V results (not now shown). It also reproduced features seen in the climatology summary of two years of data from a rooftop multi-spectral microwave radiometer in the center of Manhattan (Melecio-Vazquez et al. 2018). The seasonal mean observed θ-profiles (Fig. 18) show expected: (i) near neutral PBLs, except for the summer subsidence inversion from about 700-1200 m, (ii) higher winter variabilities, (iii) deep summer and winter daytime super-adiabatic SBLs, and  (iv) summer nighttime shallow adiabatic SBL. They also show an unexpected winter night shallow super-adiabatic SBL, probably is due to strong Q A from space heating. This result opposes the consensus since Bornstein (1968) that without synoptic influences, the entire nocturnal UBL is near adiabatic. No previous UBL study, however, had included measurements down to rooftop levels. MV21 also produced a nighttime super-adiabatic SBL, and not the near-adiabatic layer of the observations, perhaps due to an overestimated Q A . It thus is the first urban NWP model to simulate a nighttime above-rooftop super-adiabatic SBL.
The field of urban modeling has rapidly advanced so that it is now able to correctly simulate two extreme urban situations: its impacts on thunderstorms and conditions in narrow urbanized valleys. The urbanized WRF modeling study by Tomasi et al. (2019) was for Bolzano Italy, located at an altitude of 2562 m in a 2-3 km wide, deep valley. The PBL TKE scheme of Nakanishi and Niino (2009) was used with a 300 m horizontal grid spacing. Outputs were used to drive the CALPUFF dispersion model and results were compared to observations from the Bolzano Tracer Experiment (BTEX). A challenging urban weather feature to predict is when urban areas initiate or intensity a convective thunderstorm. The climatological observations of Dou et al (2015) found that an UHI of 1.5 K was necessary for this in Beijing; with lesser values existing storms bifurcated around the city. Dou et al. (2020) reproduced bifurcation for a Beijing case that used the BEP urbanization in WRF. The good results of these two studies show the progress in urban NWP modeling, some of which is based on SBL and PBL TKE results from SSZ.
Urban parameterizations requiring further development in the (BEP + BEM) + formulation, as shown by the results from the NYC study of MV21, include application of the CT formulation to all sources of Q A , not only cooling towers. The resulting additional latent heat should reduce its still over-predicted UHI values. The new variable C d formulation also requires adjustment as it over-corrects, turning large over-predicted V-values into smaller under-predicted ones. As the original formulation was based on wind tunnel simulations with uniform buildings, new simulations with the variable height building arrays typical of real cities should produce better results.
A summary report on the state of the art of urban climate research and the resiliency of cities under climate change stresses also included required future research directions (González et al. 2021). It included the following grand challenges for improved urban modeling: (i) a capability for simulation of 3-D urban PBL turbulence, (ii) representation of the nexus between urbanization, water, and energy, (iii) incorporation of human processes, (iv) prediction of urban weather and climate extremes, (v) focus on cities in coastal and complex terrain areas, and (vi) harnessing of big data. SSZ has made contributions that will enable realization of many of these goals, especially (i), (iii), and (iv), as discussed throughout.

Conclusion
This paper first summarized many scientific achievements of Sergej S. Zilitinkevich (SSZ) on the theory of surface boundary layer (SBL) and planetary boundary layer (PBL) structure. It then discussed how these theoretical advances contributed to improvements in research and operational NWP models for weather prediction, air quality (AQ), and climate change.
These include his discovery and description of his: (i) Formula for stably stratified PBL depth.
(ii) Correction to the rate equation for convectively mixed layer depth, as well as the heat and mass transfer laws for geophysical turbulent flows. (iii) Length scale of rotational stratification turbulent mixing in stably stratified PBLs.
(iv) Conceptual model of new types of atmospheric boundary layers, the: (a) conventionally neutral PBL with the stable stratification background typical of the free atmosphere that are several times thinner than truly neutral PBLs and (b) long-lived stable PBLs typical over the ocean and in winter at high latitudes. (v) Weak turbulence regime, typical of the free atmosphere, that determines the turbulent transport of energy and momentum, as well as passive scalar diffusion. (vi) 2-D horizontal TKE that forms after traditional 3-D traditional TKE is extinguished at the critical Ri. (vii) Non-local turbulent transport for by coherent structures within convective flows. (viii) Order out of chaos paradigm overturning the traditional chaos out of order paradigm.
Incorporation of these results has (and will continue to) improve the formulation and parametrizations of PBL processes in NWP, climate, and AQ models worldwide. These include use of his new.
(i) Formulations for stable layer stability functions for heat and momentum, which include both conventional 3-D TKE and the new more stable 2-D horizontal TKE, via use of the data in Figs. 2, 3, and 4. (ii) Formulations for long-lived stable marine and polar boundary layers, via the data in Fig. 1. (iii) Non-local convective mass flux schemes for unstable PBLs that capture the new "order out of chaos" insights, like that used to produce Fig. 13a. (iv) New z o and d o formulations for heat and momentum that include stability effects.
(v) Level 3 TKE closures that include a prognostic θ 2 equation and diagnostic equa-tions for horizontal turbulent fluxes. (vi) θ 2 values from (iv) for a complete heat flux formulation, as θ 2 is its counter-gradient component.
(vii) 3-D TKE formulations that include all advective and diffusion terms to correctly simulate all regions of Fig. 9. (viii) Algorithms to determine daytime PBL heights not dependent on mean quantities, but on TKE or vertical heat flux. (ix) Urban z o , d o , C d , and building TKE-production formulations that depend, not only on grid-average building heights, but also on their variations within a grid cell (x) "Personal environments" via use of crowdsourcing data within urban CFD models.
In summary, the contributions of Sergej Zilitinkevich constitute a paradigm shift in our understanding and modeling of PBL turbulence and structure, as he addressed fundamental issues such as the origin and transport properties of coherent motions, effects of buoyancy on turbulent transport, and the maintenance of turbulence in strongly stable stratification. This has produced a better understanding of atmospheric dynamics, global pollution, and climate change, allowing for future improvements in NWP, climate, and AQ models.