Toward Turbulence-Resolving Weather and Climate Simulation

(1)

Toward Turbulence-resolving

Weather and Climate Simulation

(2)

(3)

Toward Turbulence-resolving

Weather and Climate Simulation

Proefschrift

ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus prof. ir. K.C.A.M. Luyben, voorzitter van het College voor Promoties,

in het openbaar te verdedigen op dinsdag 21 april 2015 om 12:30 uur

door

Jerˆome SCHALKWIJK

ingenieur in de technische natuurkunde geboren te IJsselstein

(4)

Dit proefschrift is goedgekeurd door de promotoren: Prof. dr. H.J.J. Jonker

Prof. dr. A.P. Siebesma

Samenstelling promotiecommissie:

Rector Magnificus, voorzitter

Prof. dr. H.J.J. Jonker, promotor Prof. dr. A.P. Siebesma, promotor∗ Onafhankelijke leden:

Prof. dr. ir. H.W.J Russchenberg, Civiele Techniek en Geowetenschappen, TU Delft Prof. dr. C. Sch¨ar, ETH Z¨urich, Zwitserland

Prof. dr. B. Stevens, MPI-M, Hamburg, Duitsland

Dr. A.C.M. Beljaars, ECMWF, Reading, Verenigd Koninkrijk Prof. dr. G. van der Steenhoven, KNMI, De Bilt

Prof. dr. ir. B.J. Boersma, Technische Universiteit Delft, reservelid

∗Tevens verbonden aan het KNMI

Bij dit onderzoek is gebruik gemaakt van de supercomputer-faciliteiten van GENCI, geopereerd door CEA in het TGCC, Frankrijk, toegekend door PRACE. Daarnaast is gebruik gemaakt van de supercomputer-faciliteiten van SURF-SARA, met financi¨ele ondersteuning van NWO.

ISBN: 978-94-6295-138-9

(5)

Samenvatting

Op Weg naar Turbulentie-oplossende

Weer en Klimaat Simulatie

De vorige eeuw was getuige van de opkomst van numerieke weer- en klimaat-modellen. Weersvoorspelling, voorheen een empirische procedure gebaseerd op de persoonlijke ervaring van de meteoroloog, ontwikkelde zich tot een weten-schappelijke methode om de toekomstige toestand van de atmosfeer te voor-spellen op basis van de huidige toestand en tendens. Tegelijkertijd werden klimaatmodellen ontwikkeld om inzicht te krijgen in het gedrag van ons kli-maat over tijdspannes van meer dan honderden jaren. Als gevolg van de enorme rekenkracht die vereist is om deze weer- of klimaatmodellen te kunnen draaien, ontstond een nauwe verwikkeling tussen atmosferische en computer-wetenschap. Door de combinatie van toenemende computerkracht, beschik-baarheid van globale observaties als invoer en de onophoudelijk doorontwikke-ling van weer- en klimaatmodellen, kan de wereldwijde atmosferische stroming nu worden berekend met een ongelooflijke nauwkeurigheid.

De resterende achilleshiel van weer- en klimaatmodellen blijft echter de rep-resentatie van fenomenen op de schaal van turbulentie, zoals turbulent trans-port, wolken en neerslag. Deze fenomenen kunnen (nog) niet opgelost wor-den door hewor-dendaagse weer- en klimaatmodellen (daar hebben deze modellen te weinig detail voor), waardoor ze vertegenwoordigd moeten worden door parametrische schema’s. Om de kwaliteit van dit soort schema’s te testen zijn specialistische modellen ontwikkeld, zogenaamde Large-Eddy Simulations (LES), die een beperkt oppervlak betrekken waarop wolken wel realistisch kun-nen worden opgelost. In dit proefschrift worden LES modellen ingezet met het doel weer en klimaatvoorspellingen te verbeteren. In hoofdstuk twee worden LES resultaten gebruikt om een analytisch ’geraamte’ op te zetten met be-trekking tot de evenwichtsoplossingen van cumulus-bewolking (stapelwolken) voor verschillende grootschalige weersituaties. Zo’n opzet kan ons begrip ver-beteren van de manier waarop wolken reageren als de grootschalige situatie veranderd. Gezien weer- en klimaatmodellen projecties kunnen maken van de

(6)

ii

grootschalige situatie, kunnen we zo beter voorspellen hoe wolken zich zullen gedragen.

De ’tradionele’ manier om kennis te combineren van LES en van weer-en klimaatstudies, zoals hierbovweer-en beschrevweer-en, leunt op de mogelijkheid om de grootschalige situatie te onderscheiden van de kleinschalige situatie. Dit is tegenwoordig echter helemaal niet meer vanzelfsprekend. Computerkracht neemt dermate snel toe dat we binnen twintig jaar in staat zullen zijn om de ontwikkeling van de globale atmosfeer voor de komende week te berekenen op een turbulentie-oplossende manier: dus op een detailniveau dat het mogelijk maakt turbulentie en individuele (stapel)wolken te beschrijven. Traditionele weer- en klimaatmodellen zijn echter niet ontworpen voor dit detailniveau, en kunnen dit niet aan zonder grote ingrepen. LES modellen zijn wel ontworpen voor dit detailniveau, maar niet voor grote oppervlaktes. Het ziet er dus naar uit dat toekomstige weer- en klimaatmodellen kenmerken van beide specializa-ties zullen krijgen.

De rode draad van dit proefschrift is de studie in hoeverre we LES kun-nen gebruiken om het ’werkelijke’ weer of klimaat te modelleren. Door de grote hoeveelheid benodigde rekenkracht worden LES modellen typisch alleen toegepast voor ge¨ıdealiseerde casussen. In hoofdstuk drie wordt echter de mo-gelijkheid onderzocht om een nieuwe computertechniek toe te passen bij LES modellering: het gebruik van grafische co-processor chips (graphics processing unit, GPU). Het blijkt dat als we ons LES model aanpassen op het gebruik van deze chips, het model zodanig versneld kan worden dat nieuwe toepassingen, bijvoorbeeld in weersvoorspelling, mogelijk worden.

Het GPU-LES model wordt daarom uitgerust met procedures om op een realistische manier de grond-atmosfeer interactie en de grootschalige invloeden te behandelen. In hoofdstuk vier testen we het model door middel van simu-laties van het dagelijkse weer rondom Cabauw voor het gehele jaar 2012. Deze simulaties laten zien dat het LES model waarde kan toevoegen door zijn gede-tailleerde representatie van turbulentie en wolkenprocessen. Daarnaast wordt in hoofdstuk vijf deze simulatie gebruikt om inzicht te geven in de kwaliteit van metingen van turbulent transport van warmte en vocht. Deze metingen zijn van groot belang om een goede weerspiegeling van de atmosferische grenslaag te kunnen maken in weer- en klimaatmodellen.

Tenslotte wordt in hoofdstuk zes de grootste GPU-supercomputer van Eu-ropa ingezet om een ’proof of concept’ te bewerkstelligen van turbulentie-oplossende weersvoorspelling over Nederland met behulp van LES. Dit is een eerste kleine stap richting turbulentie-oplossende (dus met een detailniveau van individuele stapelwolken) weersvoorspelling vanuit een LES invalshoek. Het onderzoek in dit proefschrift kan zo een weg banen voor een nadere fusie van wolkensimulaties en globale weer- en klimaatmodellen.

(7)

Summary

Toward Turbulence-resolving

Weather and Climate Simulation

The previous century witnessed the advent of numerical weather and climate prediction. Weather forecasting, previously an empirical practice based on the meteorologist’ experience, became a scientific procedure to predict the future state of the atmosphere based on the current state and its tendency. Meanwhile, climate models were developed to gain insight in the behavior of the climate over more than hundreds of years. Both weather and climate modeling is highly computation-intensive, resulting in the entanglement of the atmospheric and the computer sciences. The combination of increasing computer power, the availability of global observations to serve as input and the continuous development of weather and climate models, now allows the calculation of the global atmospheric flow with incredible accuracy.

The remaining Achilles’ heel of weather and climate models, however, is the representation of turbulence scales, i.e. turbulent transport, clouds and pre-cipitation. These phenomena cannot yet be resolved by today’s global models, therefore they must be represented through parametric schemes. To test the quality of these schemes, specialized models (Large-Eddy Simulation, or LES) were created that cover only a limited area but can realistically simulate at-mospheric turbulence and clouds. This thesis aims to employ Large-Eddy Simulations to improve weather and climate predictions. In chapter two, LES results are used to build an analytic framework that describes the steady-state structure of the cumulus cloud layer for varying large-scale conditions. Such a framework may improve our understanding of the cloud response to changing conditions and improve cloud schemes. Combined with projections of changing large-scale conditions by weather and climate models, the framework may lead to an improved prediction of future weather and climate.

The ’traditional’ manner of combining knowledge from LES and from weather and climate studies, as described above, hinges on a distinction between large-scale and small-large-scale phenomena. Such a distinction is no longer self-evident,

(8)

iv

however. Computers are becoming so powerful that within the next decade or two, the global atmosphere can be integrated from today’s state to next week’s in a turbulence-resolving manner; that is, at a level of detail that al-lows the description of turbulence and individual cumulus clouds. Traditional global weather forecasting models, however, were not designed for this pur-pose, which requires the simulation of a range of small-scale phenomena. LES models, on the other hand, were designed for such a high level of detail but not for application on large areas. It seems likely that future weather and climate models will be a merger of the two approaches.

Central to this thesis is the applicability of LES modeling to the ’real’ weather and climate. Since LES models are computationally expensive, they are traditionally limited to idealized case studies. In chapter three, however, a new computational technique is explored for use in LES modeling: the utiliza-tion of the graphical co-processor chips (graphics processing unit, or GPU). It is found that utilization of these chips speeds up our LES model to such an extent that new applications, for instance in weather prediction, may become possible.

The GPU-LES model is therefore equipped to realistically treat the surface-atmosphere interaction and to handle large-scale influences. In chapter four, this model is tested for the simulation of realistic situations. In this chapter, we show results of a simulation that covers a full year of atmospheric evolution over a small area around Cabauw, the Netherlands. These simulations show that the model can provide added value through its detailed representation of turbulence and cloud processes. Chapter five shows that the model may also provide insight in, for instance, the quality of turbulent flux observations. The knowledge of turbulent fluxes is paramount to achieve an accurate estimation of surface emissions, for instance of carbon dioxide, heat and humidity, which is also essential for a good representation of the atmospheric boundary layer.

Then, in chapter six, using Europe’s largest accelerated supercomputer, the LES computation is distributed over a large number of graphical co-processors to provide a proof of concept of turbulence-resolving numerical weather pre-diction over the Netherlands. This is a first small step towards turbulence-resolving forecasting over an entire country from the angle of LES modeling. The research in this thesis thus paves the way towards a closer merger between cloud simulations and global numerical weather and climate models.

(9)

Introduction

1.1 Clouds in weather and climate

Clouds are omnipresent in our daily life. Practically every depiction of a landscape features a sky with some form of clouds, be it an urban or country landscape, be it a painting, photograph or movie. In artistic works, clouds help to capture the atmosphere of the moment, ranging from ’happy’ cumulus clouds to ominous thunderstorms.

The simple fact is that clouds do dominate our perception of the weather. A cloudless day of moderate temperature and wind feels completely different than a thickly overcast day with identical conditions otherwise. It is no sur-prise then, that many meteorological symbols on weather maps are essentially illustrations of cloudiness.

Partly, our interest in clouds is a symptom of the effects that clouds have on our habitat: the Earth’s surface. Besides being the direct source of rain, clouds have a significant temperature impact through radiative effects. These radiation effects are many: clouds reflect sunlight and absorb the Earth’s in-frared radiation, while also emitting inin-frared radiation themselves. Therefore, the net effect that clouds have on the surface is difficult to describe and pre-dict. For instance, clouds may warm the surface during nighttime but cool it during daytime.

Conversely, the temperature of the surface also influences cloud formation. A cold surface might cool the air which, as condensation occurs more easily in cold air, stimulates cloud formation. On the other hand, a cold surface supports only limited moisture evaporation and therefore may reduce the at-mosphere’s moisture content, such that cloud formation may diminish or cease altogether.

The form and coverage of clouds is thus part of an intricate balance between a large number of positive and negative feedback processes. This is one of the reasons that cloud formation is so hard to quantitatively predict. And as hard as it is to understand and predict the relative contribution of the multiple feedback processes relevant to clouds in the current climate, it is even harder to estimate the changes that will result from climate change (Stephens, 2005; Dufresne and Bony, 2008).

(14)

4 1.2. Numerical models of the atmosphere

Our understanding of cloud processes is hampered by our limited ability to perform experiments with the atmosphere. After observing the atmosphere several days at sea, we might, for instance, hypothesize that a warm sea stimu-lates cloud formation. We cannot, however, re-create the same conditions and cool the surface to test our hypothesis. We will have to wait wait for a day in which the sea surface is cooler but all other factors (humidity, wind, pressure, etc) are similar.

The entry of digital computers midway through the twentieth century has provided a partial relief. With the help of computers, we can now perform numerical experiments, or simulations, of clouds and the atmosphere.

1.2 Numerical models of the atmosphere

Clouds cannot be understood independently of the atmospheric flow of water. Moisture accounts for only a tiny fraction of the atmospheric mass (0.25%), but it plays a key role in defining the atmosphere’s characteristics (Stevens and Bony, 2013b). Clouds, in turn, make up only a tiny fraction of the atmospheric moisture, as most of the water in the atmosphere is found in the vapor state. In a sense, clouds are only the visible portion of the flow of atmospheric moisture, rather like the tip of an iceberg. The moisture flow must be understood to understand clouds. Although in principle, atmospheric motion is adequately described by the Navier-Stokes equations, these equations cannot be solved in a general manner, such that they provide only limited insight in the problem. In fact, atmospheric flow in general, and cloud formation in particular, is typically turbulent, meaning that it behaves chaotic and comprises motions on a wide range of scales. The flow of air over the globe is a superposition of motions that range between the millimeter and the planetary scale – a difference of more than a factor of 1010_{! Moreover, the motions at different scales interact} with each other, such that all scales must be known – presently obviously impossible – to get an exact solution. This problem is illustrated in the top panels of figure 1.1.

Nevertheless, solutions to the Navier-Stokes equations can be approxi-mated: we might calculate the rate of change of atmospheric air if we discretize atmospheric flow onto a number of points, or grid nodes. The scales of motion that can be resolved are related to the number of grid nodes, and the ’missing’ scales of motions are represented in a statistical fashion. Such representations, describing the statistical effect of the missing scales on the resolved scales, are called parameterizations.

The idea of numerical approximation of the weather was envisioned by Bjerknes (1904) and pioneered by Richardson (1922). Nevertheless, it took to the mid 1950s, when the first digital computers became available, for the first successful weather forecasts to be performed. In the 1960s, Lilly (1962) and Smagorinsky (1963) formulated the large-eddy simulation (LES) and the general circulation model (GCM), respectively, which would become two major lines of numerical atmospheric modeling. The length scales captured by these

(15)

1. Introduction 5 10-2 10-1 100 101 10 103 10 2 105 1064 107 108 DNS LES LAM GCM

Figure 1.1: Illustration of the length scales of atmospheric motion (in meters, cut off at 10−2 m). The green bars roughly indicate the traditional setting of the model approaches. The extra computing power available today can be used for grid refinement (blue) or domain expansion (orange).

numerical methodologies, along with later approaches, are illustrated with the horizontal bars in figure 1.1. The spectrum of atmospheric modeling can be roughly categorized as follows:

• GCMs calculate the circulation of airflow over the entire globe, and form the backbone of weather and climate forecasting. Atmospheric motions at smaller scales, traditionally including nearly all forms of clouds, must be parameterized.

• Limited-area models (LAMs), cover the middle-ground between LES (be-low) and GCM. LAM is often used as an umbrella term, describing mod-els ranging heavily parameterized circulation modmod-els on a limited area to LES-like models adapted to simulate large-scale phenomena. This category also comprises meso-scale models, e.g. the commonly used the Weather Research and Forecasting model (WRF, Skamarock et al., 2008). • LES models parameterize only the smallest scales of turbulence, which are assumed to behave universally. These models are able to realisti-cally simulate the dominant energy-containing turbulence scales. LES modeling is the approach adopted in this study.

• Direct numerical simulation (DNS ) models feature no flow parameter-ization. By solving all scales of turbulence, at the small scales these models remain closest to the laws of nature. As a result, however, they are applicable only to a small range of (mostly academic) phenomena. As computers become more powerful, the number of grid nodes (and time-steps) can increase and therefore a larger range of scales can be resolved,

(16)

6 1.3. The Gray Zone

and consequently a smaller range of scales needs to be parameterized. The computational evolution between the 1980s and today is visualized in figure 1.1, showing that the extra grid nodes can be applied in the current simulation domain to refine the solution, or they can be added at the boundaries to enlarge the simulation domain.

The advances in computational resources have directly led to advances in atmospheric science. For instance, the number of atmospheric phenomena that could be represented in large-eddy simulations grew as better computational resources became available. After Deardorff (1972) performed the case of a con-vective boundary layer, key examples are the first cumulus-topped boundary layer that was simulated by Sommeria (1976), experiments of the resolution-sensitive stratocumulus cloud by Deardorff (1980), the stable boundary layer studied by Mason and Derbyshire (1990) and the simulation of huge thunder-storm ensembles by Khairoutdinov and Randall (2006). Alternatively, more as an example of how computer development may inspire research in surprising ways, Heus et al. (2009) used a virtual reality environment to acquire a new perspective on the life-cycle of a cumulus cloud.

Other model approaches have evolved similarly. For instance, the ad-ditional computer power is applied in general circulation models to exten-sive radiation calculations, improved observation assimilation schemes, multi-simulation (ensemble) forecasting and, of course, increasing resolution (Sim-mons et al., 1989; European Centre for Medium-range Weather Forecasts, 2014).

In fact, computational advance has blurred the different modeling ap-proaches introduced in figure 1.1 up to a point at which significant overlap exists between them. The first steps have been taken to apply direct nu-merical simulations to atmospheric boundary layer flows (Jonker et al., 2013; Garcia and Mellado, 2014), the traditional domain of large-eddy simulations. Global models now feature resolutions and methods reminiscent of the tra-ditional limited-area models. Limited area models themselves can either be refined or expanded, which has led to adaptive models that behave like GCMs when employed to simulate areas that approach the Earth’s surface and like LES models when the resolution is refined (Moeng et al., 2007; Skamarock et al., 2008).

1.3 The Gray Zone

Even though computational advances may now allow it, the blending of mod-eling approaches is not as obvious as it may seem. This is due to the model designs, which are geared to a certain resolution that is related to the dis-tinction between the resolved the parameterized scales. In other words, the modeling approaches are traditionally developed to be confined to a certain range of resolutions. Nevertheless, since higher resolutions promise a better representation of both model input (e.g. topography and surface types) and hopefully also model output (e.g. precipitation), the urge to use higher

(17)

reso-1. Introduction 7

lutions has proved irresistible.

Therefore, general circulation and limited area models are now tending to resolutions that conflict with the original model design. Resulting problems arise as these models may start to resolve phenomena that are simultaneously parameterized, potentially leading to the ’double-counting’ of phenomena (e.g. Molinari and Dudek, 1992). Furthermore, parameterizations are often designed on the assumption of a sufficiently large statistical ensemble per grid box, which may be violated when the resolution is improved (e.g. Arakawa, 2004; Arakawa et al., 2011). Therefore, an increased resolution does not automatically lead to improved results. For instance, Williamson (1999) showed a case in which model output did not even converge as resolution was increased, if the param-eterizations were not modified. Likewise, Dirmeyer et al. (2012) showed that the increase of resolution had little effect on the quality of rainfall timing as long as the parameterizations remain unchanged.

These problems are particularly difficult to solve when the resolution is of the order of the energy- and flux-containing eddies (Wyngaard, 2004), since at that scale models start to resolve vertical convective overturning but still need its partial parameterization. In the atmosphere, this scale lies between the typical height of the boundary layer (∼ 1 km) and that of the troposphere (∼ 10 km). Modeling with resolutions in this range is therefore termed the ’terra incognita’ (Wyngaard, 2004) or the ’gray zone’. The problems associated with the gray zone are acute, since global numerical weather prediction models, now operating at resolutions of roughly 16 km (European Centre for Medium-range Weather Forecasts, 2014), are rapidly approaching it and regional models, at resolutions of the order of several kilometers (e.g. Steppeler et al., 2003; Lean et al., 2008; Seity et al., 2011), have already entered the gray zone.

Arakawa et al. (2011) identified two possible approaches to tackle the gray zone problem. One approach would be to construct a unified parameterization that is ’scale-aware’, that is, it adapts its behavior to the employed resolution. However, the task to a construct such a unified parameterization is immense. A second approach can therefore be identified, which may be investigated complementary to the first and involves the use of cloud-resolving models to act as parameterizations. Several strategies to this approach have been proposed (e.g. Grabowski, 2001; Jung and Arakawa, 2010) with a similar general idea: cloud-resolving models are run sparsely in space to provide the information on the small scales that otherwise requires parameterization. As computing power increases, the role of cloud-resolving models can gradually increase up to the point at which they cover the planet. Part of the research in this thesis can be used to assess whether Large-Eddy Simulations are suitable as cloud-resolving – even turbulence-resolving – models in such a framework.

1.4 Computer facilities

Besides the parameterization challenges related to resolution increase, the us-age of extra computing power itself is also seldom effortless. None of the

(18)

nu-8 1.5. This thesis

merical experiments mentioned above were performed using a standard desk-top computer. As science tends to concern itself with the boundaries of our understanding, so do numerical studies tend to use the maximally available computational resources. As a result, numerical scientists had to adapt to changing methods with which these supercomputers were instructed.

Early digital computers were instructed using punch cards, which were labor-intensive to create or edit. When magnetic tape storage became avail-able, programming became easier and programming languages became more versatile. As supercomputer facilities developed that consisted of an increasing number of processing units, or cores, programmers had to adopt parallel algo-rithms and protocols. One of the recent changes in (super)computers, affecting both hardware and software, is the introduction of the graphics processing unit (GPU) for general computing purposes.

The GPU originates from the consumer computer market. Towards the end of the twentieth century, the human-computer interface became increasingly graphical in nature. Computers were equipped with operating systems that allowed visual insight and manipulation of data. Since the central processing unit (CPU) became increasingly burdened with rendering graphics, a relatively monotonous task, it became economical to employ computers with a dedicated, auxiliary processing unit dedicated to that task: the GPU. This unburdened the CPU to provide a more rapid computing experience.

As graphics rendering became more complex and detailed, the GPU was developed to cope with the increasing demands. In fact, due to its special-ized architecture, the GPU started to bypass the CPU in terms of computing power. When this was appreciated, the general-purpose GPU was developed: a GPU that could be instructed for general computations (i.e. other than graph-ics rendering). General computations, in this context, must be taken with a grain of salt: the strength of the GPU lies in its dedication to the monotonous task of matrix and vector computations associated with graphics rendering, making it suitable mostly for tasks that are mathematically relatively similar. More specifically, tasks that are well-suited for the GPU are (1) computation-intensive, (2) parallel and (3) favor throughput over latency (i.e. the speed of the bulk of computations is more important than the speed of any indi-vidual computation) (Owens et al., 2008). For applications which meet these characteristics, the GPU is highly cost- and energy-efficient. Therefore, since these characteristics describe a broad range of problems in high-performance computing, a large number of modern supercomputer systems now also feature GPU units.

1.5 This thesis

It turns out that atmospheric computations in general, and large-eddy simula-tions in particular, are ideally suited for the GPU. This was keenly recognized by Griffith (2010), who took on the audacious step to adapt the Dutch Atmo-spheric Large-Eddy Simulation (DALES; developed by Nieuwstadt and Brost,

(19)

1. Introduction 9

1986; Cuijpers and Duynkerke, 1993; Heus et al., 2010, among many others) to utilize the GPU. Part of the aim of this research is to continue the develop-ment of the resulting model and make it suitable for modern GPU-accelerated supercomputer systems.

Although the development of a numerical model is not a goal in itself, the GPU offers advantages that directly translate to new possibilities for numerical experiments, for instance:

• On just a single GPU, situated in a desktop computer, numerical ex-periments can be performed that were ’state-of-the-art’ only a few years back, and still worthwhile today. These simulations can be performed in a graphical manner that is not altogether different from gaming: sim-ulations are visualized and can be edited graphically before and during the simulation. The direct visual feedback has obvious advantages in feasibility studies and case set-up, but is also promising for teaching purposes.

• Secondly, single-GPU simulations can evade the relatively slow computer-to-computer communication associated with traditional parallel simula-tions, while still enjoying the many-core GPU architecture. The result is that small-scale simulations (up to, say, 2563 _{grid nodes) can be} per-formed much quicker than they could by using parallel CPU systems, which are more adapted to simulations with many grid nodes. This im-plies that such simulations can be numerically integrated over a much longer time than otherwise feasible.

• Finally the traditional, and most obvious advantage: the possibility to simulate more grid nodes, and thus an increased range of turbulent scales. By coupling a large number of GPUs, extraordinarily large simulations can be performed without requiring an extraordinary amount of com-puter cores.

The possibilities outlined here recur in the work in this thesis:

Chapter 2 scans a phase-space of large-scale conditions by performing a large array of simulations,

Chapter 3 shows the possibilities of interactive simulation, Chapters 4 and 5 are the result of a very long simulation, and Chapter 6 describes a very large-area simulation.

1.6 Outline

An overarching theme of LES-based research as introduced above is the im-provement of our capability to provide weather and climate projections. The traditional method is to study, and formulate simplified frameworks for, e.g. cloud formation using LES. An attempt to formulate a simple but realistic analytic framework to non-precipitating cumulus convection is presented in

(20)

10 1.6. Outline

chapter 2. Here the aim is to develop a modeling approach for clear and cloudy convection that is as simple as possible while retaining the essential behavior, such that the equilibrium response of the cloudy boundary layer to varying climatological forcings can be examined. The analytic results are guided and verified by a large ensemble of simulations, made possible by an easy visual case set-up and quick results provided by the GPU version of the LES. In fact, the framework developed in this chapter can be used to set up a generic case of steady-state cumulus convection under arbitrarily chosen (within specific limits) large-scale conditions. This framework might aid further research in establishing the changes in cumulus cloud layers in a changing climate.

Chapter 3 further explains the advantages to atmospheric GPU computing and introduces GALES: a GPU-resident Atmospheric Large-Eddy Simulation. GALES, initiated by Griffith (2010), was brought level to its CPU counter-part and equipped with an interactive, graphical user interface that allows researches to interactively influence running simulations. In this chapter, the benefits of GPU acceleration in particular, and interactive graphical simula-tions in general, to numerical atmospheric research are discussed.

Keeping figure 1.1 in mind, it becomes apparent that the advantages pro-vided by GALES, as outlined in the previous section, also allow for a more di-rect applications of LES in weather and climate prediction. Such applications may include 1) small-scale ’nowcasts’: high-resolution, short-term forecasts over a selected area, 2) ’super-parameterizations’ in which LES models take on the role of the parameterizations in a GCM (as in, e.g. Grabowski, 2001) or even 3) directly employing LES to provide turbulence-resolving limited-area weather prediction.

Chapter 4 may prove a relevant exploratory exercise of the first two pos-sibilities. It presents a year-long large-eddy simulation run that simulates the weather over Cabauw, the Netherlands, by taking input from a KNMI limited area model in the framework of the KNMI Parameterization Testbed (Neggers et al., 2012). This exercise is unique in that it represents, to the author’s knowledge, the first continuous year-long Large-Eddy Simulation run. The results of this run thus comprise a wide variety of weather conditions and in-cludes numerous transitions, at a location at which also detailed observations are available. Therefore, the resulting dataset can provide an opportunity to study the statistics of a multitude of weather phenomena in the LES. The re-sults are compared with observations to validate and understand the behavior of the LES model.

Chapter 5 follows up on this chapter and presents a detailed analysis of the turbulent fluxes in the year-long run. The year-long data-set provides insight in which time-scales are dominant for turbulent transport, and the implications on attempts to quantify turbulent transport from single-point observations.

Although furthest into the future still, a proof of concept of the principle of using LES models to provide limited-area forecasts is given in chapter 6, in which a 100m-resolution large-eddy simulation of the weather of the Nether-lands is presented. This chapter shows a hint of how numerical weather and climate forecasting may look ’beyond’ the gray zone, when turbulence, clouds

(21)

1. Introduction 11

and convection are explicitly resolved. It is also attempted to provide a broader perspective by speculating on possible approaches to future weather prediction.

(22)

(23)

CHAPTER 2

Simple Solutions to Steady-state

Cumulus Regimes in the Convective

Boundary Layer

A modeling framework is developed which extends the mixed layer model to steady state cumulus convection. The aim is to consider the simplest model which retains the essential behavior of cumulus-capped layers. The presented framework allows for the evaluation of stationary states depen-dent on external parameters. These states are completely independepen-dent of the initial conditions, and therefore represent an asymptote which might help deepen our understanding of the dependence of the cloudy boundary layer on external forcings. Formulating separate equations for the lifting condensation level and the mixed layer height, the dry and wet energet-ics can be distinguished. Comparing the dry and wet buoyancy effects, regimes are identified which can support steady-state cumulus clouds, and regimes which cannot. The dominant mechanisms are identified which govern the creation and eventual depth of the cloud layer. Model pre-dictions are tested by comparison with a large number of independent Large-Eddy Simulations, for varying surface and large-scale conditions, and are found to be in good agreement.

2.1 Introduction

Boundary layer clouds play an important role in both the dynamical and ra-diative properties of the boundary layer, controlling to an important extent the height, the efficiency of vertical transport and the opacity of the boundary layer. However, these clouds are also notoriously hard to model due to the high resolution needed to resolve boundary layer turbulence. Stratocumulus cloud

(24)

14 2.2. Bulk modeling

decks require high vertical resolution to resolve cloud-top entrainment, and properties of cumulus clouds are also sensitive to the horizontal grid resolution due to their inhomogeneity.

Understanding of cumulus clouds has increased significantly over the past decades, in part due to Large Eddy Simulations (LES, e.g. Deardorff, 1970), which allowed many numerical studies to be performed on detailed cloud pro-cesses. Especially understanding the cloud-environment interactions (Paluch, 1979; Reuter and Yau, 1987; Siebesma and Cuijpers, 1995; Heus and Jonker, 2008) has been a popular topic over many years, enhancing our understanding of cumulus clouds and improving parameterizations. However, the behavior of the cumulus-capped boundary layer as a whole and the interplay with large scale tendencies and forcings is still poorly understood. This is further em-phasized by studies on the wide spread of cloud-climate feedbacks (Dufresne and Bony, 2008; Bony et al., 2006), revealing the large uncertainties associated with the feedback behavior of low clouds.

For this reason, this work attempts a modeling approach as simple as pos-sible while retaining the essential behavior. Starting with the mixed-layer model for clear convective situations, building on the work of Tennekes (1973); Lilly (1968) and, more recently, Vil`a-Guerau de Arellano (2004); van Driel and Jonker (2011), the model is further elaborated to include a cloud layer, work-ing in the line of Stevens (2006), Neggers et al. (2006), Stevens (2007), Bellon and Stevens (2012) and Bellon and Stevens (2013). The goal of this exercise is to set up a framework in which one can study the response of cumuliform clouds to large-scale atmospheric forcings and surface properties, minimizing complications by simplifying the system to its essence. Note that in doing this, a somewhat different approach is chosen than, for example, Bretherton and Park (2008), Nuijens and Stevens (2012), who propose more realistic mod-els that focus on the dynamical temporal response of the cloud layer, whereas we focus on the stationary solutions.

The focus on stationary solutions allows us to formulate analytical solutions for the regimes where cumuliform clouds are supported. These solutions can be solved analytically or numerically, but do not require model integration in time. The hope is that such solutions further promote our understanding of the limiting factors on cumulus regimes. Moreover, these solutions can be used to quickly set up numerical experiments of steady-state cumulus convection in varying environmental conditions and forcings, which might aid future steady-state cumulus studies.

2.2 Bulk modeling

Consider a conserved scalar ψ. The evolution of the ensemble average ψ is described by the following conservation equation (e.g. Stevens, 2006):

∂ψ ∂t + w ∂ψ ∂z = − ∂φψ ∂z + Sψ (2.1)

(25)

2. Simple Solutions to Steady-state Cumulus Regimes in the CBL 15 ψm ψ₊ ψf(z) ∆ψ ψ0 h a) ψ(z) z Mixed layer ψm ψ₊ ψf(z) ∆ψ ψ0 b) h L ψ(z) Cloud layer Mixed layer

Figure 2.1: Assumed profile of a conserved scalar ψ for the case of a) a clear boundary layer and b) cumulus convection. In this framework, the ’jump’ is considered the difference between the value above the boundary layer (+) and the mixed-layer value. Note that the profile for cumulus convection collapses to the clear boundary layer when L > h.

assuming incompressibility. In this equation, φψ denotes the vertical turbu-lent flux of scalar ψ, i.e. φψ = w′ψ′, where primes indicate deviations from the ensemble mean. Sources, including radiative cooling and horizontal advec-tion, are contained in the source term Sψ, which can depend on height. The ensemble-mean vertical velocity w(z) is associated with large-scale motions.

In this manuscript we will focus on the steady-state solutions, allowing us to reduce equation (2.1) to:

w∂ψ ∂z = −

∂φψ

∂z + Sψ (2.2)

Now assume the boundary layer to be described by profiles as depicted in figure 2.1. The profiles are assumed to be well mixed throughout the mixed layer. The mixed layer reaches up to level L – the lifting condensation level (LCL) – which marks the start of the cloud layer. The boundary layer top is denoted as h, and we assume turbulent fluxes to vanish at this level. The cloud layer is thus included in h. In the event no clouds form, the profile is well-mixed all the way up to the boundary layer height h. Above h, ψ is described by the free tropospheric state ψf.

Boundary conditions are provided by the free troposphere above the bound-ary layer height and the surface value ψ0 in an infinitesimally thin surface

(26)

layer. We integrate over this profile from z = 0 to z = h+_{, where we define} h+_{= lim}

ǫ↓0h + ǫ, just above boundary layer top, to incorporate the disconti-nuity at the top into the integration.

Using the mixed layer assumption that ψ(z) = ψm for 0 < z < min(L, h), the integration yields:

w(h) ψ+− ψm − Z h+ L ψ(z) − ψm ∂w(z) ∂z dz = φψ,0+ Z h+ 0 Sψdz (2.3)

where the turbulent flux at h+ _{vanishes. Subscripts are used as shorthand} notation for the evaluation point, i.e. φψ,0 represents the surface flux and ψ+ is short for ψ(h+_{) = ψ}f_(h).

The effects of an eventual cloud layer gradient enter the equation through the the second term on the left hand side. In a recent study, Bellon and Stevens (2012) attempted to capture these cloud layer effects in the parameter γ, which can be reformulated in the current framework as follows:

γ = 1 − 1 w(h) Z h+ L ψ(z) − ψm ψ+− ψm ∂w(z) ∂z dz (2.4)

which captures the effect of the interaction between the large-scale vertical velocity w and the cloud layer gradient in a single dimensionless quantity. Using (2.4) one can rewrite equation (2.3) as

γw(h) ψ+− ψm= φψ,0+ Z h+

0

Sψdz (2.5)

Bellon and Stevens (2012) simplified (2.5) by assuming γ to have a constant value of 0.8. To simplify even further, in this study we ignore the cloud layer gradient, which boils down to setting γ = 1; the sensitivity of the results to γ will be evaluated in section 2.42.4.3.

2.2.1 Boundary layer height

In order to determine the steady-state boundary layer height h, we consider the boundary layer energetics. We follow the approach of Stevens (2007) and use what is sometimes called the ’liquid water virtual potential temperature’ (Grenier and Bretherton, 2001).

θvl= θl+ ǫIθqt (2.6)

with ǫI = R_Rv_d − 1. In the absence of liquid water, θvl is equal to the virtual potential temperature. As the temperature variation in the second term is typically negligible, hereafter the liquid water virtual potential temperature is approximated as a linear combination of θland qt, i.e. eǫI ≡ ǫIθ ≈ constant.

(27)

2. Simple Solutions to Steady-state Cumulus Regimes in the CBL 17

The result is that, to good approximation, θvl is conserved, implying that equation (2.1) also holds for θvl. In the mixed-layer, ∂θvl/∂z = 0 and equation (2.2) for θvl reduces to:

∂φθvl(z)

∂z = Sθvl(z) for z < min(L, h) (2.7)

In steady state, the flux divergence is required to balance the sources. If φθvl,0

is known, equation (2.7) can be integrated to determine φθvl(z) throughout

the mixed layer.

This allows us to define the ’dry thermal reach’ η, which describes the height that convection would reach were there no effects from condensation (i.e. ”dry” in the sense that no liquid water effects are taken into account). We define η as the height at which:

φθvl(η

−_{) = −aφ}

θvl(0) (2.8)

with a a constant fraction, representing the ’entrainment efficiency’. Equa-tion (2.8) describes one of the the fundamental properties of buoyancy-driven boundary layers (Ball, 1960; Betts, 1973): the energy provided by the surface buoyancy flux dissipates during the upward motion in such a way that we find the top of the cloudless boundary layer at h = η in equation (2.8). Hence, equation (2.8) represents the classical mixed-layer closure in case no clouds form, with η = h. In case clouds form, η becomes a maximum bound on the mixed-layer height, as we must have L < η for clouds to form.

The value of η can be solved from combining equations (2.7–2.8), which provides an implicit equation for η:

−φθvl,0(1 + a) =

Z η−

0

Sθvl(z)dz (2.9)

Note that it is η which provides the criterion for cumulus cloud formation: clouds will form if the lifting condensation level is found below η, i.e. for L ≤ η, whereas the boundary layer remains cloudless when L > η. The lifting condensation level L is determined by the mixed-layer values of θl,mand qt,m, as well as the surface pressure ps,

L = f(θl,m, qt,m, ps) (2.10) If L > η, the steady-state solution describing the boundary layer state is found from equations (2.5) and (2.9) with ψ = θvland h = η.

2.2.2 Cloud layer

In general, equation (2.5) conveys that θl,m and qt,m will be a function of the boundary layer height h. From equation (2.10), then so is the lifting condensation level: L = L(h). Hence, a cloudless solution is possible if L(η) > η, i.e. the thermodynamic properties of a boundary layer with height h = η are such that L > η and no clouds will form.

(28)

In case L(η) ≤ η, cloud formation will occur and an alternative equation for h is needed, since h 6= η. In typical bulk models, the system is then closed using the assumption that L = η. Much success has been achieved using this assumption, and we will start by using this closure as well.

It is worthwhile to investigate the premises for this closure. The equation for η is based on the height dry thermals reach, whereas L is the lifting conden-sation level following from θl,mand qt,m. This difference between η and L was already recognized by Betts and Ridgway (1989). Why would dry thermals be depleted of their energy exactly at the lifting condensation level?

Consider a situation where L(η) < η, i.e. η − L(η) > 0. Thermals will now reach the lifting condensation level, and form clouds from that level on-ward. The extra energy from the release of latent heat allows the boundary layer depth h to increase beyond the ’dry thermal reach’ η. As a result the thermodynamic properties of the layer are modified, which in turn lead to a modified lifting condensation level L′ = L(h) and, in the most general case, to a modified η′.

Two possibilities can be distinguished. If η′− L′ < η − L, the lifting con-densation level approached the dry thermal reach. This allows fewer thermals to reach this level, thereby limiting the total latent heat release. In the ex-treme event that L′> η′, no thermals will reach L′ and no latent heat release follows. This is a self-correcting mechanism with the classical closure as limit:

L = η (2.11)

As L is a function of h, equation (2.11) can be regarded an implicit equation for the steady-state boundary layer height h.

The second possibility, however, is that η′− L′> η − L. This would allow more thermals to reach lifting condensation level, increasing latent heat release. This increase would result in a further increased h, which in turn results in η′′− L′′ > η′− L′, etc. This situation is unstable. It becomes apparent that to exclude such instability, we must require that

∂η − L

∂h < 0 (2.12)

at some height h for a stable steady-state boundary layer with that height to be possible.

In conclusion, we need equation (2.12) as a premise to close the system using L = η. Moreover, this closure is a steady-state result of the adaptation of the boundary layer to the thermodynamic effect of latent heat release. For that reason, the use of this closure in a time-dependent framework implicitly assumes that the timescales of this feedback are much smaller than the other timescales in the system.

2.2.3 Free atmosphere

In the current framework we also assume the free troposphere to be in steady state, which entails that equation (2.2) also applies to ψf(z) for all z > h.

(29)

Because the turbulent fluxes are zero for z > h, and because the equation must hold for any value of h, one arrives at the following relation between the source Sψ(z), the subsidence profile, and the free tropospheric profile ψ

f (z) w(z)∂ψ f (z) ∂z = Sψ(z) (2.13)

The sources (advection, radiation, etc.) balance the effect of subsidence through-out the atmosphere, independent of h.

2.2.4 Solutions

Given ψ ∈ {θl, qt}, the steady-state solution can now be summarized as follows. Equation (2.5) can be solved for the mixed-layer values of θl and qt. In our aim to simplify the model to its essence, we start by setting γ = 1. This allows a smooth transition from cloudless to cloudy regimes. The result is:

θl = θlf(h) − φθl,0+ Rh+ 0 Sθldz w(h) (2.14) qt = qft(h) − φqt,0+ Rh+ 0 Sqtdz w(h) (2.15)

Note that we have dropped the overbars as well as the subscript m on {θl, qt} for readability. Unless explicitly mentioned otherwise, hereafter θl and qt de-note θl,m and qt,m, respectively.

The dry thermal reach η acts as a bound to determine whether a cloud layer appears. It can be solved from the implicit equation (2.9), repeated here in order to compactly present all model equations:

−φθvl,0(1 + a) =

Z η−

0

Sθvl(z)dz (2.16)

Provided ∂L/∂h > 0, the boundary layer height h can be solved from equation (2.11)

h = (

η if L(η) > η

L−1(η) if L(η) ≤ η (2.17)

In other words, if clouds form, the boundary layer height is found by solving L(h) = η for h. In the absence of clouds, h = η and the solution reduces to the mixed layer equations for the clear boundary layer as formulated by Tennekes (1973) and still used in recent studies (Vil`a-Guerau de Arellano, 2004; van Driel and Jonker, 2011), when these are solved for steady state. In that case, the equation dh/dt = we_{+ w(h) = 0 requires the entrainment rate w}e _to balance the subsidence, which is directly implied in equations (2.14–2.16).

(30)

20 2.3. Model settings θf0_l θl,m θl,+ L h θl z qt qtf0 qt,m qt,+ Free atmosphere Cloud layer Mixed layer

Figure 2.2: Schematic picture of the connection between the idealized ’free atmosphere’ profile, extending from above the boundary layer all the way down to the ground (dashed lines), and the actual profile of the atmosphere, with a boundary layer up to z = h (solid line). The boundary layer is divided into a well-mixed part and a cloudy part.

2.3 Model settings

2.3.1 Case set-up

In order to describe the environment in a simplified yet generic way, we define the free atmospheric state as the profile the troposphere would have, would there be no boundary layer (extrapolating downward). Now consider the fol-lowing idealized free atmospheric profiles (dashed lines in figure 2.2)

θf_l(z) = θf 0_l + Γz (2.18)

qft(z) = q f 0

t (2.19)

The lapse rate Γ is kept constant over height and in time. For given w(z), equations (2.18–2.19) also determine the sources through equation (2.13).

The formation of the boundary layer can alter these profiles only up to the boundary layer height h such that a top boundary condition is automatically supplied. The resulting ’jumps’ ∆ψ = ψ+− ψm= ψf(h) − ψmare defined by ψmand h as follows

∆θl = θlf 0+ Γh − θl,m (2.20)

(31)

Table 2.1: Environmental parameters used for the phase-space exploration. Parameter Value

θf 0_l 290 K

qtf 0 0 kgkg−1

Γ 6 × 10−3 Km−1

ps 102900 P a

A schematic picture of this notion is drawn in figure 2.2. The values of the pa-rameters Γ, θf 0_l and qtf 0are provided in table 2.1. All Large-Eddy Simulations performed in this study start with initial profiles as described by equations (2.18–2.19), i.e. θl(z) = θ_lf(z) and qt(z) = qft(z) at t = 0.

By prescribing the free atmospheric conditions and the subsidence profile (as we will), we have defined the sources through equation (2.13), i.e. Sqt = 0

and

Sθl= wΓ (2.22)

2.3.2 Parameter evaluation

The entrainment efficiency a sets the maximum height dry thermals reach and thus controls the boundary layer height of the cloudless boundary layer. Therefore, we diagnose a from all cloudless Large-Eddy Simulations (15 in total) presented in section 2.42.4.1–2.42.4.2. The LES model is described in appendix A. The boundary layer height the simulations reach is diagnosed, after which a is chosen such that equation (2.16) best captures the trend.

On the basis of those data, we have found a to be best described by a value of 0.4. For this diagnosis, we defined η as the location of the maximum gradient in the virtual potential temperature profile, following Sullivan et al. (1998). While this value for a might seem large, it represents the best zero order representation of the boundary layer profile in cloudless cases. It does not represent the minimum of the buoyancy flux, which will always be smaller in magnitude. Note that a is the only control parameter in this zero order model.

A short discussion on the numerical value of a is provided in appendix B.

2.3.3 Phase-space explorations

In this study we will investigate the behavior of the boundary layer for different external forcings. In particular, we will consider the following cases.

Case I: In section 2.42.4.1, the model is studied in case of a constant diver-gence D, i.e. w(z) = −Dz. To begin we explore the simple case of fixed surface fluxes, in which we study the influence of their magnitude and composition.

(32)

22 2.4. Results

Case II: In section 2.42.4.2, we explore the effect of coupled surface fluxes, varying the sea surface temperature and D. The effect of the simplifying assumption to set γ = 1 will be considered in section 2.42.4.3. We will investigate the influence of γ and the physical mechanism it captures, reconsidering the results of section 2.42.4.2.

Case III: In section 2.42.4.4 the subsidence is chosen constant with height. Among others, this allows us to test the relevance of the assumption that ∂L/∂h > 0, since ∂L/∂h = 0 in this case.

Finally, in section 2.5 the dominant mechanisms will be identified and the behavior in different forcing regimes will be compared.

2.4 Results

2.4.1 Case I: Fixed fluxes and divergence

To begin we explore the simple case of fixed surface fluxes in case of constant divergence D, i.e. w(z) = −Dz. In this case, equations (2.13–2.16) can be evaluated as: θl(h) = θf 0_l +φθl,0 Dh + 1 2Γh (2.23) qt(h) = qf 0t + φqt,0 Dh (2.24)

with the dry thermal reach η:

η =

r

2φθvl,0(1 + a)

DΓ (2.25)

Note that in this case, ∂η/∂h = 0. Now, the dependence of the lifting condensation level on h can be investigated by expanding:

∂L ∂h = ∂L ∂qt ∂qt ∂h + ∂L ∂θl ∂θl ∂h (2.26)

where we know that ∂L/∂qt < 0 and ∂L/∂θl > 0 from thermodynamical arguments. Therefore, substituting equations (2.23–2.24) and using that h > η, equation (2.26) can be used to show that:

∂L

∂h > 0 (2.27)

for φθvl,0 > 0. This justifies the use of equation (2.17), which together with

equations (2.23–2.25) describes the steady-state solution of {θl, qt, h, L} if the parameters {φθl,0, φqt,0, θ

f 0 l , q

f 0

(33)

−0.04 −0.020 0 0.02 0.04 10 20 30 40 Bowen ratio β [-] S u rf a ce b u o y a n cy fl u x B0 [W /m 2] a) ζ <50m ζ ≥50m −0.04 −0.020 0 0.02 0.04 500 1000 1500 Bowen ratio β [-] C lo u d d ep th ζ [m ] b) −0.04 −0.02 0 0.02 0.04 500 1000 1500 2000 Bowen ratio β [-] B o u n d a ry la y er h ei g h t h [m ] c) B0= 10W/m2 B₀= 15W/m2 B0= 20W/m2 B0= 25W/m2 293K < T0<299K −0.04 −0.02 0 0.02 0.04 500 1000 1500 Bowen ratio β [-] η [m ] d) −0.04 −0.02 0 0.02 0.04 290 292 294 296 298 Bowen ratio β [-] θl [K ] e) −0.04 −0.024 0 0.02 0.04 6 8 10 12 14 Bowen ratio β [-] qt [g /k g ] f)

Figure 2.3: Steady state solutions for parameters as in table 2.1 in a phase-space of surface buoyancy flux B0 and Bowen ratio β. The shaded contours in a) illustrate the qualitative behavior of the model prediction for cloud layer depth ζ = h − L, the black line representing the separation between cloudy and cloudless regimes. The overlaid symbols each depict the state of a LES simulation after 12 days of simulated time, shape and fill color representing the LES’ cloud layer depth. Time series of the the cloud layer depth for the case of B0 = 20W m−2 are shown in b). Model results for cloud layer depth, mixed-layer height, mixed-layer θl and qtare outlined in more detail through lines for selected values of B0 in panels c)–f ), while the symbols depict LES results. The grey lines illustrate the path which is traversed if fluxes are set interactively as a function of the sea surface temperature.

(34)

24 2.4. Results

To illustrate these solutions in a phase-space of external parameters, we describe the surface fluxes in terms of the surface buoyancy flux B0≡ ρcpφθvl,0

(W m−2) and the Bowen ratio β

β = cp Lv

φθl,0

φqt,0

(2.28)

with cpthe specific heat capacity of water and Lv the latent heat of vaporiza-tion. The fluxes are prescribed and have a constant value over time.

The choice of surface fluxes allows a simple interpretation throughout the phase-space. The surface buoyancy flux B0describes how much kinetic energy is brought into the boundary layer from the surface, and thereby governs η. Note that equation (2.25) is a function of the surface buoyancy flux only; so for a given surface buoyancy flux, η is independent of the Bowen ratio.

The Bowen ratio determines how the surface energy is divided over latent and sensible heat. Therefore, it governs the values of θl and qt and thus the critical lifting condensation level L(η). A lower Bowen ratio corresponds to a moister and cooler boundary layer, while η remains constant. This obviously corresponds to a situation more favorable to cloud formation.

The results are shown in figure 2.3. Panel a illustrates the model solution by depicting the cloud layer depth ζ ≡ h − L using shaded contours. The so-lutions are compared with a number of Large-Eddy Simulations (LES). Each LES starts with initial profiles equal to the idealized free atmospheric state described by equations (2.18–2.19), thus without an initial boundary layer. Environmental parameters are provided in table 2.1, and we prescribed a di-vergence of D = 7 × 106_s₋₁_{. The LES simulations are further described in the} appendix.

The simulations are performed for 12 days to allow the simulations to reach a steady-state. The symbols in figure 2.3a each represent a separate LES simulation, circles denote cases with ζ < 50m and squares denote cases with a steady ζ ≥ 50m; their color indicates the cloud layer depth. As we study the phase-space in the close vicinity of the cloudless/cloudy boundary, LES simulations in the ’cloudless’ regime will still feature some especially moist updrafts which reach their condensation point. Therefore, we chose ζ = 50m as the boundary for LES simulations, as simulations with less than that will typically have a significantly lower liquid water path. The transition from clear to cloudy is illustrated more quantitatively in panel c.

Simulations which either ran into the top of the domain or did not reach steady state are depicted with a plus (+) and are not considered in the panels c–f. This only occurs for very thick cloud layers and is likely to be an artifact of the LES set-up: re-performing random samples of such cases with a larger (vertical) domain and integrating over longer times, these simulations do reach well-defined steady-states.

A simulation is considered to have reached a steady-state if the 6hr-mean altitude ztop, defined as the lowest altitude where qt(ztop) ≤ 1 × 10−5kg kg−1, remains steady over at least 24 hours.

(35)

Each LES simulation ends in its own state, closely corresponding to the model predictions, while each started from the same initial conditions. Fig-ure 2.3b shows the time series of cloud depth for selected simulations with B0= 20 W m−2. This panel shows how cases with large Bowen ratios remain cloudless, while those with decreasing Bowen ratios develop deeper cloud lay-ers, while each simulation started without clouds altogether. The simulation with β = −0.02 ran into the domain top, but the time series does suggest converging behavior.

Given the surface buoyancy flux, the model predicts a clear boundary layer at equal height for all large enough Bowen ratios. Clouds form as the Bowen ratio decreases, increasing the boundary layer height. The mixed-layer values of θl and qt change accordingly, changing L up to the point that L = η. Considering that η in equation (2.25) is a function of surface buoyancy flux only, the surprising result is that the steady-state lifting condensation level is determined by the surface buoyancy flux, and independent of Bowen ratio.

This behavior is compared with LES results in more detail in panels c–d. Although the cloud layer depth is systematically underestimated by the model, the trend is similar to LES. LES simulations with equal surface buoyancy flux show remarkably little variation in the mixed-layer height, which is diagnosed by considering the average cloud-base height or – in cloudless cases – the maximum gradient in the virtual potential temperature.

The underestimation of cloud layer depth is most likely to result from neglecting the cloud layer gradient. A further discussion on this subject is provided in section 2.42.4.3.

Panels e–f show the mixed-layer values of the temperature and humidity. Especially in humidity, the effect of cloud formation is readily observed. For low Bowen ratios, the moistening of the mixed layer is inhibited by the fact that clouds form, which ventilates the excess moisture into the cloud layer (Neggers et al., 2006). Cloud formation simultaneously increases the entrainment of dry and warm air into the boundary layer, explaining the warming for low Bowen ratios in panel e.

2.4.2 Case II: Fixed sea-surface temperature

One could argue that the approach of assuming constant surface fluxes is un-realistic, even in the stationary limit. For one, such an approach neglects the feedback which occurs as the result of the difference between the boundary layer state and the surface properties. An interesting approach, therefore, is to allow this feedback to occur by introducing limited surface-atmosphere interaction as follows

φθl,0 = V (θl,0− θl) (2.29)

φqt,0 = V (qt,0− qt) (2.30)

where θl,0 is the liquid water potential temperature related to the sea surface temperature T0via surface pressure psand qt,0the saturation humidity at T0.

(36)

26 2.4. Results

In this approach, we have simplified the transfer coefficient V = CD|U| as a constant to avoid feedbacks in the velocity U and to allow for non-zero surface fluxes while remaining in the idealized case of no mean winds. We set V to the value of 1 cm/s, which would correspond to the transfer coefficient of a case with typical values of wind speed U = 5 m/s and CD= 2 × 10−3.

When this surface parameterization is introduced into equations (2.14– 2.15), the steady-state solutions become

θl(h) = Dhθ f 0 l + V θl,0+1₂ΓDh2 Dh + V (2.31) qt(h) = Dhq f 0 t + V qt,0 Dh + V (2.32)

revealing a explicit combination of surface and top conditions acting on the boundary layer. The surface-buoyancy flux reaches a steady-state solution as follows φθvl,0(h) = DhV Dh + V θvl,0− θvlf 0− 1 2Γh (2.33) where θvl,0= θl,0+ ˜ǫIqt,0 represents the sea surface virtual potential temper-ature.

The largest increase in complexity with regard to the previous section ap-pears in the solution for the layer height η. The solution for the mixed-layer height, equation (2.25), remains valid:

η(h) = r

2φθvl,0(h)(1 + a)

DΓ (2.34)

but since the magnitude of the surface buoyancy flux is dependent on h, so is η. In case of cloudless convection, though, equation (2.34) is an implicit equation for η since then h = η.

Phase-space

The coupling between cloud layer and mixed layer has become further inter-twined through air-surface interaction. Whereas the boundary layer height h already depended on the mixed-layer dynamics in case I, η was independent of the boundary layer height until now. The result is that none of the model variables {θl, qt, L, h} can be regarded as independent of the other.

Nevertheless, when steady-state has been reached, fluxes are constant and as such, the system is equivalent to the system with constant surface fluxes. Essentially, we have lost a degree of freedom in the set up of the cases, since the surface fluxes are now governed by a single parameter – the sea surface temperature T0 – instead of two. Varying the sea surface temperature, all other parameters held constant, thus amounts to following a trajectory through the {β, B0} phase-space in figure 2.3. This path is visualized as a gray line in all panels of figure 2.3, representing the steady-state solution in case of an interactive surface flux found by varying T0 from 293K to 299K. These

(37)

2 4 6 8 10 292 294 296 298 300 S ea su rf a ce te m p er a tu re [K ] D [10−6_s−1_] a) ζ <50m ζ ≥50m 2 4 6 8 10 12 0 200 400 600 800 1000 Time [days] C lo u d d ep th ζ [m ] b) 2 4 6 8 10 0 500 1000 1500 D [10−6_s−1_] C lo u d d ep th ζ [m ] c) T0= 293K T0= 295K T0= 297K T0= 299K 2 4 6 8 10 0 500 1000 1500 2000 D [10−6_s−1_] η [m ] d) 2 4 6 8 10 290 292 294 296 298 D [10−6_s−1_] θl [K ] e) 2 4 6 8 10 8 10 12 14 16 D [10−6_s−1_] qt [g /k g ] f)

Figure 2.4: Steady state solutions for parameters as in table 2.1. Shaded contours in a) depict the model for cloud layer depth ζ = h − L, and overlaid symbols depict LES simulations as before. Time series of the the cloud layer depth for T0 = 295K are shown in b). Model results for cloud layer depth, mixed-layer height, mixed-layer θl and qt are outlined in more detail in line plots for selected values of T0 in panels c)–f ).

(38)

28 2.4. Results

solutions are compared to four Large-Eddy Simulations shown by star-shaped markers. The line starts at β = −0.03 and B0 = 5W m−2 at T0 = 293K. At this point evaporation drives the boundary layer to such an extent that the entrainment of warm air causes the boundary layer to become warmer than the sea surface, resulting in a negative sensible heat flux. Similar behavior was observed by Nuijens and Stevens (2012). As T0 increases, β and B0 increase to about 0.03 and 30W m−2, respectively.

For a fair comparison, the LES sea surface feedback was simplified to be-have as equations (2.29–2.30), with the mixed-layer values replaced by the first model level. The simulations are not found exactly on the gray line as they are free to develop surface fluxes which might differ from the model solutions. They are very close, however, demonstrating how an interactive flux mecha-nism moves through the surface flux phase-space. Note that the surface fluxes traverse the phase-space in a line nearly parallel to the boundary separating cloudy and cloudless regimes.

Figure 2.3a also illustrates the reason we chose to study the regime of rather low Bowen ratios, −0.04 ≤ β ≤ 0.04, since this is where the regime change in varying T0 occurs. For comparison, the BOMEX intercomparison case (Siebesma et al., 2003) has a Bowen ratio of 0.06 (but has non-comparable free-tropospheric conditions).

Since we have lost a degree of freedom in the surface fluxes, we can now more easily study the behavior as a function of divergence D, as well. We have constructed a phase-space of T0 and D in figure 2.4 with equivalent panels as figure 2.3. The line through the surface flux phase-space of figure 2.3 is now the vertical line in figure 2.4a at D = 7 × 10−6s−1. The choice to scan a phase-space of sea surface temperature and subsidence allows us to perform all LES simulations, like before, starting from the exact same initial condi-tions. The sea surface temperature and divergence are varied; the cooling by radiation and large-scale horizontal advection varies such that the divergence is always balanced (equation 2.13). This implies that the sink term increases in amplitude with D, i.e. the cooling is reduced if D decreases and vice-versa. Even though the added feedbacks increase the complexity in the dynam-ics of the system, the general behavior can be understood using the notions from previous sections. An increased sea surface temperature results in an in-crease in B0 and β. Given constant free atmospheric conditions, an increased sea surface temperature therefore leads to a warmer, moister, and therefore higher boundary layer. A warmer sea favors cloud formation as the increased buoyancy allows thermals to reach the lifting condensation level more easily.

At the same time, the dry thermal reach typically increases as divergence decreases, corresponding to intuition. Cumulus clouds, therefore, form more easily in regions of lower divergence, as this allows the mixed-layer to deepen and reach the lifting condensation level.

Typical parameter values for the trade wind areas in the current climate (Bony et al., 2006; Zhang et al., 2009) correspond to the upper left corner of figure 2.4a). It should be no surprise that this region is nested firmly in the cumulus regime.

Toward Turbulence-Resolving Weather and Climate Simulation

Toward Turbulence-resolving

Weather and Climate Simulation

Toward Turbulence-resolving

Weather and Climate Simulation

Proefschrift

Samenvatting

Op Weg naar Turbulentie-oplossende

Weer en Klimaat Simulatie

Summary

Toward Turbulence-resolving

Weather and Climate Simulation

Contents

CHAPTER 1

Introduction

1.1

Clouds in weather and climate

1.2

Numerical models of the atmosphere

1.3

The Gray Zone

1.4

Computer facilities

1.5

This thesis

1.6

Outline

CHAPTER 2

Simple Solutions to Steady-state

Cumulus Regimes in the Convective

Boundary Layer

2.1

Introduction

2.2

Bulk modeling

2.2.1

Boundary layer height

2.2.2

Cloud layer

2.2.3

Free atmosphere

2.2.4

Solutions

2.3

Model settings

2.3.1

Case set-up

2.3.2

Parameter evaluation

2.3.3

Phase-space explorations

2.4

Results

2.4.1

Case I: Fixed fluxes and divergence

2.4.2

Case II: Fixed sea-surface temperature