We consider the closure problem of representing the higher-order moments (HOMs) in terms of lower- order moments, a central feature in turbulence modeling based on the Reynolds-averaged Navier–Stokes (RANS) approach. Our focus is on models suited for the description of asymmetric, nonlocal, and semiorganized turbulence in the dry atmospheric convective boundary layer (CBL). We establish a multivariate probability density function (PDF) describ- ing populations of plumes that are embedded in a sea of weaker randomly spaced eddies, and apply an assumed delta-PDF approximation. The main content of this approach consists of capturing the bulk properties of the PDF. We solve the clo- sure problem analytically for all relevant HOMs involving velocity components and temperature and establish a hierarchy of new non-Gaussian turbulence closure models of different content and complexity ranging from analytical to semianalyti- cal. All HOMs in the hierarchy have a universal and simple functional form. They refine the widely used Millionshchikov closure hypothesis and generalize the famous quadratic skewness–kurtosis relationship to higher order. We examine the performance of the new closures by comparison with measurement, LES, and DNS data and derive empirical constants for semianalytical models, which are best for practical applications. We show that the new models have a good skill in predict- ing the HOMs for atmospheric CBL. Our closures can be implemented in second-, third-, and fourth-order RANS turbu- lence closure models of bi-, tri-, and four-variate levels of complexity. Finally, several possible generalizations of our approach are discussed.