Adaptive Cartesian meshes for atmospheric single-column models: a study using Basilisk 18-02-16

It is well known that the representation of certain atmospheric conditions in climate and weather models can still suffer from the limited grid resolution that is facilitated by modern-day computer systems. Herein we study a simple one-dimensional analogy to those models by using a single-column model description of the atmosphere. The model employs an adaptive Cartesian mesh that applies a high-resolution mesh only when and where it is required. The so-called adaptive-grid model is described, and we report our findings obtained for tests to evaluate the representation of the atmospheric boundary layer, based on the first two GEWEX ABL Study (GABLS) inter-comparison cases. The analysis shows that the adaptive-grid algorithm is indeed able to dynamically coarsen and refine the numerical grid whilst maintaining an accurate solution. This is an interesting result as in reality, transitional dynamics (e.g. due to the diurnal cycle or due to changing synoptic conditions) are the rule rather than the exception.

1 Introduction

Single-column models (SCMs) are often used as the building blocks for global (or general) circulation models (GCMs). As such, many of the lessons learned from SCM development can be inherited by GCMs and hence the evaluations of SCMs receive considerable attention by the geoscientific model development community (Baas et al., 2017; Bosveld et al., 2014; Neggers et al., 2012). In this work, we present an SCM that employs an adaptive Cartesian mesh that can drastically reduce the computational costs of such models, especially when pushing the model's resolution. The philosophy is inspired by recently obtained results on the evolution of atmospheric turbulence in a daytime boundary layer using three-dimensional (3-D) adaptive grids. As promising results were obtained for turbulence-resolving techniques such as direct numerical simulations and large-eddy simulation (LESs), herein we explore whether similar advancements can be made with more practically oriented techniques for the numerical modelling of the atmosphere. As such, the present model uses Reynolds-averaged Navier–Stokes (RANS) techniques to parameterize the vertical mixing processes due to turbulence (Reynolds, 1895), as is typical in weather and climate models.

The discussion of limited grid resolution is present in many studies of SCMs and GCMs. A prominent example is the nocturnal cumulus-cloud case (Wyant et al., 2007): whereas a high-resolution mesh is required to capture the processes at the cloud interface, a coarser resolution may be used for the time when the sun has risen and the cloud has been dissolved. More generally speaking, virtually all of the atmospheric dynamics that require a relatively high-resolution grid for their representation in numerical models are localized in both space and in time. The issue is made more difficult to tackle by the fact that their spatio-temporal localization is typically not known a priori (e.g. the height and strength of a future inversion layer). Therefore, the pre-tuned and static-type grids that most operational GCMs use (virtually all) are not flexible enough to capture all dynamical regimes accurately or efficiently. This also puts a large strain on the closures used for the sub-grid-scale processes. In order to mitigate this challenge, GCMs that employ a so-called adaptive grid have been explored in the literature. Here, the grid resolution adaptively varies in both space and time, focussing the computational resources on where and when they are most necessary. Most notably, the innovative work of Jablonowski (2004), Jablonowski et al. (2009), and St-Cyr et al. (2008) report on the use of both Cartesian and non-conforming three-dimensional adaptive grids and clearly demonstrate the potential of grid adaptivity for GCMs. Inspired by their work, we follow a 1-D SCM approach and aim to add to their findings, using different grid-adaptative formulations and solver strategies. Since SCMs do not resolve large-scale atmospheric circulations, the analysis herein focusses on the representation of the atmospheric boundary layer (ABL).

Over the years, the computational resources that are available to run computer models have increased considerably (Schaller, 1997). This has facilitated GCMs to increase their models' spatial resolution, enabling us to resolve the most demanding processes with increased grid resolutions. However, it is important to realize that the (spatial and temporal) fraction of the domain that benefits most from an increasing maximum resolution necessarily decreases as a separation of the modelled spatial scales increases (Popinet, 2011). This is because the physical processes that warrant a higher-resolution mesh are virtually never space filling. For example, the formation phase of tropical cyclones is localized in both space and time and is characterized by internal dynamics that evolve during the formation process. By definition, with an increasing scale separation, only an adaptive-grid approach is able to reflect the effective (or so-called fractal) dimension of the physical system in the scaling of the computational costs (Popinet, 2011; Van Hooft et al., 2018). This is an aspect where the present adaptive-grid approach differs from, for example, a dynamic-grid approach (Dunbar et al., 2008), which employs a fixed number of grid cells that needs to be predefined by the user. This work employs a similar method for grid adaptation as presented in the work of Van Hooft et al. (2018) on 3-D turbulence-resolving simulations of the ABL. As such, this work is also based on the adaptive-grid toolbox and built-in solvers provided by the “Basilisk” code (http://basilisk.fr, last access: 25 November 2018).

We test our model with the well-established cases defined for the first two Global Energy and Water cycle EXchanges (GEWEX) ABL Study (GABLS) inter-comparison projects for SCMs. As part of the GEWEX modelling and prediction panel, GABLS was initiated in 2001 to improve our understanding of the atmospheric boundary layer processes and their representation in models. Based on observations during field campaigns, a variety of model cases has been designed and studied using both LES and SCMs with a large set of models using traditional static-grid structures. An overview of the results and their interpretation for the first three inter-comparison cases are presented in the work of Holtslag et al. (2013). Here, we will test the present adaptive-grid SCM based on the first two inter-comparison cases, referred to as GABLS1 and GABLS2. These cases were designed to study the model representation of the stable boundary layer and the diurnal cycle, respectively. Their scenarios prescribe idealized atmospheric conditions and lack the complete set of physical processes and interactions encountered in reality. At this stage within our research, the authors consider this aspect to be an advantage, as the present SCM model does not have a complete set of parameterizations for all processes that are typically found in the operational models (see, e.g., Grell et al., 2005; Slingo, 1987).

This paper is organized as follows. The present SCM is discussed in more detail in Sect. 2. Based on the results from a simplified flow problem, Sect. 3 starts with an analysis of the numerical methods used and the grid-adaptation strategy. Model results for ABL-focussed cases that are based on the first two GABLS inter-comparison scenarios are also presented in Sect. 3. Finally, a discussion and conclusions are presented in Sect. 4.

2 Model overview

As we focus on the merits of grid adaptivity in this study on SCMs and not on the state-of-the-art closures for the vertical transport phenomena, we have opted to employ simple and well-known descriptions for the turbulent transport processes. More specifically, the present model uses a stability-dependent, first-order, local K-diffusivity closure as presented in the work of Louis et al. (1982) and Holtslag and Boville (1993). For the surface-flux parameterizations we again follow the formulations in the work of Holtslag and Boville (1993). However, to improve the representation of mixing under stable conditions, an alteration is made to the formulation of the so-called stability correction function under stably stratified conditions. Based on the work of England and McNider (1995), we use a so-called short-tail mixing function. The closures used for the turbulent transport are summarized next. The upward surface fluxes (F) of the horizontal velocity components (u, v), the potential temperature (θ), and specific humidity q are evaluated as

$$\begin{array}{}\text{(1a)}& {\displaystyle }& {\displaystyle }{F}_u=-C_M{U}_{\mathrm{1}}{u}_{\mathrm{1}},\text{(1b)}& {\displaystyle }& {\displaystyle }{F}_v=-C_M{U}_{\mathrm{1}}{v}_{\mathrm{1}},\text{(1c)}& {\displaystyle }& {\displaystyle }{F}_{\mathit{\theta }}=-C_H{U}_{\mathrm{1}}\left({\mathit{\theta }}_{\mathrm{1}}-{\mathit{\theta }}_{\mathrm{0}}\right),\text{(1d)}& {\displaystyle }& {\displaystyle }{F}_q=-C_H{U}_{\mathrm{1}}\left(q_{\mathrm{1}}-q_{\mathrm{0}}\right),\end{array}$$

where U is the wind-speed magnitude and indices 0 and 1 refer to the values at the surface and the first model level, respectively. The surface transport coefficients are

$$\begin{array}{}\text{(2a)}& {\displaystyle }& {\displaystyle }{C}_M=C_{\mathrm{N}}{f}_{\mathrm{s},M}\left({\mathit{\text{Ri}}}_{\mathrm{b}}\right),\text{(2b)}& {\displaystyle }& {\displaystyle }{C}_H=C_{\mathrm{N}}{f}_{\mathrm{s},H}\left({\mathit{\text{Ri}}}_{\mathrm{b}}\right),\end{array}$$

with Ri_b the surface bulk Richardson number, which is defined as

$$\begin{array}{}\text{(3)}& {\displaystyle }{\mathit{\text{Ri}}}_{\mathrm{b}}={\displaystyle \frac{g}{{\mathit{\theta }}_{\mathrm{v},\mathrm{ref}}}}{\displaystyle \frac{z_{\mathrm{1}}\left({\mathit{\theta }}_{\mathrm{v},\mathrm{1}}-{\mathit{\theta }}_{\mathrm{v},\mathrm{0}}\right)}{U_{\mathrm{1}}^{\mathrm{2}}}},\end{array}$$

where g is the acceleration due to gravity, θ_v is the virtual potential temperature, and θ_v,ref is a reference temperature whose value is taken as a scenario-specific constant. Equation (3) assumes that θ_v is related to the buoyancy (b) (Boussinesq, 1897) via g and θ_v, ref according to $b=g/{\mathit{\theta }}_{\mathrm{v},\mathrm{ref}}({\mathit{\theta }}_v-{\mathit{\theta }}_{\mathrm{v},\mathrm{ref}})$ . The virtual potential temperature is related to the potential temperature (θ) and specific humidity (q) according to

$$\begin{array}{}\text{(4)}& {\displaystyle }{\mathit{\theta }}_{\mathrm{v}}=\mathit{\theta }\left(\mathrm{1}-\left(\mathrm{1}-{\displaystyle \frac{R_{\mathrm{v}}}{R_{\mathrm{d}}}}\right)q\right),\end{array}$$

with $R_{\mathrm{v}}/R_{\mathrm{d}}=\mathrm{1.61}$ the ratio of the gas constants for water vapour and dry air (Emauel, 1994; Heus et al., 2010). The so-called neutral exchange coefficient (C_N) is calculated using

$$\begin{array}{}\text{(5)}& {\displaystyle }{C}_{\mathrm{N}}={\displaystyle \frac{k^{\mathrm{2}}}{\ln {\left((z_{\mathrm{1}}+z_{\mathrm{0},M})/z_{\mathrm{0},M}\right)}^{\mathrm{2}}}},\end{array}$$

with k=0.4 the von Kármán constant, z₁ the height of the lowest model level, and z_0,M the roughness length for momentum. For the cases studied in this work, the roughness length for heat is assumed to be identical to z_0,M. The stability correction functions for the surface transport of momentum and heat (f_s, M, f_s, H) are

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{f}_{\mathrm{s},M}\left({\mathit{\text{Ri}}}_{\mathrm{b}}\right)\\ \text{(6a)}& {\displaystyle }& {\displaystyle }=\left\{\begin{array}{ll}\mathrm{0},& {\mathit{\text{Ri}}}_{\mathrm{b}}\ge \mathrm{0.2},\\ {\left(\mathrm{1}-\frac{{\mathit{\text{Ri}}}_{\mathrm{b}}}{\mathrm{0.2}}\right)}^{\mathrm{2}},& \mathrm{0}\le {\mathit{\text{Ri}}}_{\mathrm{b}}<\mathrm{0.2},\\ \mathrm{1}-\frac{\mathrm{10}{\mathit{\text{Ri}}}_{\mathrm{b}}}{\mathrm{1}+\mathrm{75}{C}_{\mathrm{N}}\sqrt{\left((z_{\mathrm{1}}+z_{\mathrm{0},M})/z_{\mathrm{0},M}\right)\Vert {\mathit{\text{Ri}}}_{\mathrm{b}}\Vert }},& {\mathit{\text{Ri}}}_{\mathrm{b}}<\mathrm{0},\end{array}\right.{\displaystyle }& {\displaystyle }{f}_{\mathrm{s},H}\left({\mathit{\text{Ri}}}_{\mathrm{b}}\right)\\ \text{(6b)}& {\displaystyle }& {\displaystyle }=\left\{\begin{array}{ll}{f}_{\mathrm{s},M}\left({\mathit{\text{Ri}}}_{\mathrm{b}}\right),& {\mathit{\text{Ri}}}_{\mathrm{b}}\ge \mathrm{0},\\ \mathrm{1}-\frac{\mathrm{15}{\mathit{\text{Ri}}}_{\mathrm{b}}}{\mathrm{1}+\mathrm{75}{C}_{\mathrm{N}}\sqrt{\left((z_{\mathrm{1}}+z_{\mathrm{0},M})/z_{\mathrm{0},M}\right)\Vert {\mathit{\text{Ri}}}_{\mathrm{b}}\Vert }},& {\mathit{\text{Ri}}}_{\mathrm{b}}<\mathrm{0},\end{array}\right.\end{array}$$

which conclude the description of the surface fluxes. The vertical flux ($\stackrel{\mathrm{‾}}{w^{\prime }{a}^{\prime }}$) of a dummy variable a due to turbulence within the boundary layer is based on a local diffusion scheme and is expressed as

$$\begin{array}{}\text{(7)}& {\displaystyle }\stackrel{\mathrm{‾}}{w^{\prime }{a}^{\prime }}=-K{\displaystyle \frac{\partial a}{\partial z}},\end{array}$$

where K is the so-called eddy diffusivity

$$\begin{array}{}\text{(8)}& {\displaystyle }K=l^{\mathrm{2}}Sf\left(\mathit{\text{Ri}}\right).\end{array}$$

l represents an effective mixing length,

$$\begin{array}{}\text{(9)}& {\displaystyle }l=\min \left(kz,l_{\mathrm{bl}}\right),\end{array}$$

with l_bl being the Blackadar length scale; we use, l_bl=70 m (Holtslag and Boville, 1993). S is the local wind-shear magnitude,

$$\begin{array}{}\text{(10)}& {\displaystyle }S=\sqrt{{\left({\displaystyle \frac{\partial u}{\partial z}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial v}{\partial z}}\right)}^{\mathrm{2}}},\end{array}$$

and f(Ri) is the stability correction function for the vertical flux,

$$\begin{array}{}\text{(11)}& {\displaystyle }f\left(\mathit{\text{Ri}}\right)=\left\{\begin{array}{ll}\mathrm{0},& \mathit{\text{Ri}}\ge \mathrm{0.2},\\ {\left(\mathrm{1}-\frac{\mathit{\text{Ri}}}{\mathrm{0.2}}\right)}^{\mathrm{2}},& \mathrm{0}\le \mathit{\text{Ri}}<\mathrm{0.2},\\ \sqrt{\mathrm{1}-\mathrm{18}\mathit{\text{Ri}}},& \mathit{\text{Ri}}<\mathrm{0},\end{array}\right.\end{array}$$

i.e. based on the gradient Richardson number,

$$\begin{array}{}\text{(12)}& {\displaystyle }\mathit{\text{Ri}}={\displaystyle \frac{g}{{\mathit{\theta }}_{\mathrm{v},\mathrm{ref}}}}{\displaystyle \frac{\partial {\mathit{\theta }}_v/\partial z}{S^{\mathrm{2}}}}.\end{array}$$

The authors of this work realize that there have been considerable advancements in the representation of mixing under unstable conditions in the past decades, e.g non-local mixing (Holtslag and Boville, 1993) and closures based on turbulent kinetic energy (Lenderink and Holtslag, 2004; Mellor and Yamada, 1982). Therefore, we would like to note that such schemes are compatible with the adaptive-grid approach, and they could be readily employed to improve the physical descriptions in the present model. From an implementations' perspective, those schemes would not require any grid-adaptation-specific considerations when using the Basilisk code.

For time integration, we recognize a reaction–diffusion-type equation describing the evolution of the horizontal wind components and scalar fields such as the virtual potential temperature and specific humidity (q). For a variable field s(z,t), we write,

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\partial s}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(K{\displaystyle \frac{\partial }{\partial z}}s\right)+r,\end{array}$$

where r is a source term and K is the diffusion coefficient; c.f. Eq. (8). Using a mixed implicit–explicit first-order-accurate time discretization for the diffusive term and an explicit time integration for the source term (r), with time step Δt separating the solution sⁿ and sⁿ⁺¹, this can be written

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{s^{n+\mathrm{1}}-s^n}{\mathrm{\Delta }t}}={\displaystyle \frac{\partial }{\partial z}}\left(K^n{\displaystyle \frac{\partial }{\partial z}}{s}^{n+\mathrm{1}}\right)+r^n.\end{array}$$

Rearranging the terms we write

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\partial }{\partial z}}\left(K^n{\displaystyle \frac{\partial }{\partial z}}{s}^{n+\mathrm{1}}\right)-{\displaystyle \frac{s^{n+\mathrm{1}}}{\mathrm{\Delta }t}}=-{\displaystyle \frac{s^n}{\mathrm{\Delta }t}}-r^n\end{array}$$

to obtain a Poisson–Helmholtz equation for sⁿ⁺¹, using the eddy diffusivity calculated from the solution sⁿ (Kⁿ). Eq. (15) is solved using a multigrid strategy, employing a finite-volume-type second-order-accurate spatial discretization (Popinet, 2017a, b). The source term r in Eq. (13) is defined using different formulations for the various scalar fields in our model. For θ and q, the source term r concerns the tendency in the lowest grid level due to the surface fluxes (F, see Eqs. 1, $r_{F_s}$) and the effect of large-scale synoptic divergence (r_w) according to the vertical velocity w (i.e. prescribed for the GABLS2 case). We write, for a dummy variable s,

$$\begin{array}{}\text{(16)}& {\displaystyle }{r}_{w,s}=-w{\displaystyle \frac{\partial s}{\partial z}}.\end{array}$$

For the horizontal velocity components (u, v) the corresponding source terms (i.e. $r_{F_s}$ and r_w) are also taken into account and supplemented with the additional source term ${\mathit{r}}_{{\mathbf{\nabla }}_{\mathit{h}}P,f}$, which concerns the horizontal pressure-gradient-forcing vector (i.e. $-\frac{\mathrm{1}}{\mathit{\rho }}{\mathbf{\nabla }}_{\mathit{h}}P$, for air with a density ρ) and the Coriolis-force term according to the local Coriolis parameter f. For the horizontal velocity vector $\mathit{u}=\mathit{\{}u,v,\mathrm{0}\mathit{\}}$ we write

$$\begin{array}{}\text{(17)}& {\displaystyle }{\mathit{r}}_{{\mathbf{\nabla }}_{\mathrm{h}}P,f}={\displaystyle \frac{-{\mathbf{\nabla }}_{\mathrm{h}}P}{\mathit{\rho }}}-f\left(\widehat{\mathit{k}}\times \mathit{u}\right),\end{array}$$

where “×” represents the cross product operator and $\widehat{\mathit{k}}=\mathit{\{}\mathrm{0},\mathrm{0},\mathrm{1}\mathit{\}}$ the unit vector in the vertical direction. In this work we adopt the commonly used strategy to introduce a velocity vector that is known as the geostrophic wind (U_geo), according to

$$\begin{array}{}\text{(18)}& {\displaystyle }{\mathit{U}}_{\mathrm{geo}}={\displaystyle \frac{\widehat{\mathit{k}}}{\mathit{\rho }f}}\times {\mathbf{\nabla }}_{\mathrm{h}}P.\end{array}$$

The most prominent feature of the SCM presented in this work is that it adaptively coarsens and refines the grid resolution based on the evolution of the solution itself. As mentioned in the introduction, the associated grid-adaptation algorithm is the same as described in Van Hooft et al. (2018). Here, we only briefly discuss the general concept.

Apart from the imperfect representation of the physical aspects of a system in numerical models, additional errors naturally arise due to the spatial and temporal discretization. In general, a finer resolution corresponds to a more accurate solution and a simulation result is considered to be “converged” when the numerically obtained solution and the statistics of interest do not crucially depend on the chosen resolution. The aim of the grid-adaptation algorithm is to dynamically coarsen and refine the mesh so that the errors due to the spatial discretization remain within limited bounds and are uniformly distributed in both space and time. For our adaptive approach this requires (1) an algorithm that evaluates a local estimate of the discretization error in the representation of selected solution fields (χ_a for a field “a”) and (2) a corresponding error threshold (ζ_a) which determines if a grid cell's resolution is either too coarse (i.e. χ_a>ζ_a), too fine (i.e. ${\mathit{\chi }}_a<\mathrm{2}{\mathit{\zeta }}_a/\mathrm{3}$), or just fine. Grid adaptation can then be carried out accordingly and the solution values of new grid cells can be found using interpolation techniques. A cell is refined when the estimated error for at least one selected solution field exceeds its refinement criterion and a cell is coarsened when it is considered to be “too fine” for all selected solution fields. The “error estimator” (χ) is based on a so-called multi-resolution analysis that is formally linked to wavelet thresholding. The algorithm aims to estimate the magnitude of higher-order contributions in the spatial variability of the solution that are not captured by the solver's numerical schemes. Consistent with the second-order spatial accuracy of the solver's numerical schemes (Popinet, 2017b), the algorithm employs a second-order accurate wavelet-based error estimate. In practice, grid refinement will typically occur at the locations where the solution is highly “curved”, whereas those areas where the solution fields vary more “linearly” in space are prone to coarsening. The error threshold, or so-called refinement criterion ζ, is defined by the user. Noting that similar to the pre-tuning of the fixed-in-time grids as is common in most SCMs, the balance between accuracy and the required computational effort remains at the discretion of the model's user.

For the cases in this work that focus on the ABL (i.e. in Sects. 3.2 and 3.3), the dynamics are governed by the wind ($\mathit{U}=(u,v)$) and the virtual potential temperature (θ_v); hence we base the refinement and coarsening of the grid on a second-order-accurate estimated error associated with the representation of these discretized fields. Based on trial and error, we set the corresponding refinement thresholds

$$\begin{array}{}\text{(19)}& {\displaystyle }{\mathit{\zeta }}_{u,v}& {\displaystyle }=\mathrm{0.25}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}},\text{(20)}& {\displaystyle }{\mathit{\zeta }}_{{\mathit{\theta }}_{\mathrm{v}}}& {\displaystyle }=\mathrm{0.5}\mathrm{K}\end{array}$$

for both of the horizontal wind components and virtual potential temperature, respectively. These values are the result of a choice by the authors that aims to strike an arbitrary balance between the accuracy of the solution and the computational effort required to run the model. Note that a similar (arbitrary) balance needs also to be found when static grids are employed. For a simple flow set-up, Sect. 3.1 presents an example convergence study to show the effects of using different refinement criteria on the accuracy of the obtained solutions.

Table 1 The exact formulation of the methods are described at the online locations of the definition files for the different cases presented in this paper.

The wall-clock times are evaluated using a single core (processor model: Intel i7-6700 HQ).

Grid adaptation is carried out each time step. The tree-based anisotropic-grid structure in Basilisk facilitates a convenient basis for the multi-resolution analysis and the subsequent refinement and coarsening of cells at integer levels of refinement. This entails that the spatial resolution can vary by factors of 2 (Popinet, 2011). For the adaptive-grid runs presented in this paper, the time spent in the actual grid assessment and adaptation routines is less than 5 % of the total wall-clock time (see Table 1).

Apart from the Ekman-spiral case in Sect. 3.1, the physical time step in the ABL-focussed cases is adaptively varied between 2 and 15 s based on the convergence properties of the aforementioned iterative solver. Note that these values are rather small compared to existing GCMs that often employ higher-order-accurate time-integration schemes. Additionally, the correlation of spatial and temporal scales warrants a smaller time step, since the present model employs a higher maximum vertical resolution compared to that of an operational GCM. The solver's second-order spatial accuracy is validated and the performance is accessed for a simple flow set-up in Sect. 3.1. For the exact details of the model set-ups for the cases presented in this paper, the reader is referred to the case-definition files (in legible formatting). Links are provided to their online locations in Table 1.

3 Results

3.1 The laminar Ekman spiral and grid adaptation

Before we focus our attention on cases that concern the ABL, this section discusses the philosophy of the grid-adaptation strategy used based on the analysis of a one-dimensional (1-D) laminar Ekman-flow set-up. This simple and clean set-up enables us to quantify numerical errors explicitly and test the solver's numerical schemes. The aim of this section is to show that the grid-adaptation strategy and the accompanying refinement criteria provide a consistent and powerful framework for adaptive mesh-element-size selection. Results are presented for both an equidistant-grid and the adaptive-grid approach. The case describes a neutrally stratified fluid with a constant diffusivity for momentum (K) given by the kinematic viscosity ν and density ρ in a rotating frame of reference with respect to the Coriolis parameter f. A flow is forced by a horizontal pressure gradient (−∇_hP) according to Eq. (18) using ${\mathit{U}}_{\mathrm{geo}}=\left\{U_{\mathrm{geo}},\mathrm{0}\right\}$, over a no-slip bottom boundary (located at z_bottom=0). Assuming that the velocity components converge towards the geostrophic wind vector for z→∞ and vanish at the bottom boundary, there exists an analytical, 1-D, steady solution for the horizontal wind component profiles (u_E(z), v_E(z)), which is known as the Ekman spiral;

$$\begin{array}{}\text{(21)}& {\displaystyle }{u}_{\mathrm{E}}& {\displaystyle }=U_{\mathrm{geo}}\left(\mathrm{1}-e^{-\mathit{\gamma }z}\cos \left(\mathit{\gamma }z\right)\right),\text{(22)}& {\displaystyle }{v}_{\mathrm{E}}& {\displaystyle }=U_{\mathrm{geo}}{e}^{-\mathit{\gamma }z}\sin \left(\mathit{\gamma }z\right),\end{array}$$

with γ the so-called inverse Ekman depth; $\mathit{\gamma }=\sqrt{f/\left(\mathrm{2}\mathit{\nu }\right)}$. We choose numerical values for U_geo,γ, and f of unity in our set-up and present the results in a dimensionless framework. The solution is initialized according to the exact solution, and we set boundary conditions based on Eqs. (21) and (22). Equation (13) is solved numerically for both u and v components, on a domain with height $z_{\mathrm{top}}=\mathrm{100}{\mathit{\gamma }}^{-\mathrm{1}}$. The simulation is run until $t_{\mathrm{end}}=\mathrm{10}{f}^{-\mathrm{1}}$, using a fixed time step $\mathrm{\Delta }t=\mathrm{0.01}{f}^{-\mathrm{1}}$. The time step is chosen sufficiently small such that the numerical errors are dominated by the spatial discretization rather than by the time-integration scheme. During the simulation run, discretization errors alter the numerical solution from its exact, and analytically steady, initialization. For all runs, the diagnosed statistics regarding the numerical solutions that are presented in this section have become steady at t=t_end.

Figure 1 The locally evaluated error in the numerical solutions for u and v at t=t_end, based on the analytical solution (ϵ_u,v, see Eq. 23) for 10 runs using different equidistant mesh-element sizes. Panel (a) shows that the diagnosed errors for each run plotted against the mesh-element size used (Δ) times the inverse Ekman depth (γ, see text). Panel (b) shows, with the same colour coding as in (a), the correlation between the wavelet-based estimated error (χ) and the corresponding diagnosed error in the numerically obtained solution (ϵ). The inset (using the same axis scales) shows the results for a single run and reveals a spread of several orders of magnitude in both ϵ and χ values.

Figure 2 The locally evaluated error in the numerical solutions for u and v at t=t_end, based on the analytical solution (ϵ_u,v, see Eq. 23) for 20 runs using the adaptive-grid approach with different refinement criteria (see colour bar). Panel (a) shows that the diagnosed errors for each run plotted against the mesh-element size used (Δ). The inset (using the same axis scales) shows the results for a single run. Panel (b) shows the correlation between the wavelet-based estimated error (χ) and the corresponding diagnosed error in the numerically obtained solution (ϵ). The inset (using the same axis scales) shows the results for a single run and reveals a relatively small spread in both ϵ and χ values compared to the equidistant-grid results presented in Fig. 1b.

Figure 3 The scaling characteristics for the laminar Ekman-spiral case. Panel (a) presents the error convergence for the equidistant-grid and adaptive-grid approach. The errors (η) follow the slope of the blue dashed line that indicates the second-order accuracy of the methods. The wall-clock time for the different runs is presented in (b), showing that for both of the aforementioned approaches, the required effort scales linearly with the number of grid cells.

The spatial-convergence properties for the equidistant-grid solver are studied by iteratively decreasing the (equidistant) mesh-element sizes used (Δ) by factors of 2, and we monitor the increasing fidelity of the solution at t=t_end between the runs. Therefore, based on the analytical solution, a local error (ϵ_u,v) of the numerically obtained solution (uⁿ, vⁿ) within each grid cell is diagnosed and is defined here as

$$\begin{array}{}\text{(23)}& {\displaystyle }{\mathit{ϵ}}_a=\left|a^n-〈a_{\mathrm{E}}〉\right|,\end{array}$$

where a is a dummy variable for u and v, 〈a_E〉 is the grid-cell-averaged value of the analytical solution (a_E), and aⁿ is the value of the numerical solution within the cell. Note that aⁿ also represents the grid-cell-averaged value in our finite-volume approach. Figure 1a shows the results for all runs and compares the grid resolution used (Δ) with the error ϵ_u,v. It appears that the observed range of ϵ values is large and typically spans 10 orders of magnitude, with a lower bound defined by the “machine precision” (i.e. $\approx {\mathrm{10}}^{-\mathrm{15}}$ for double-precision floating-point numbers). This wide range can be explained by the fact that the Ekman spiral is characterized by exponentially decreasing variation with height (see Eqs. 21, 22), and hence the equidistant grid may be considered overly refined at large z. This illustrates that, for a given solver formulation, the error in the solution is not directly dictated by the mesh-element size but also depends on the local shape of the numerical solution itself. This poses a challenge for the pre-tuning of meshes applied to GCMs, where a balance needs to be found between accuracy and computational speed performance. The solution of a future model run is not known beforehand, and hence the tuning of the grid typically relies heavily on experience, empiricism, and a priori knowledge. This motivates us to apply the method of error estimation in the representation of a discretized solution field as described in Popinet (2011) and Van Hooft et al. (2018). For both velocity components, this estimated error (χ_u,v) is evaluated at the end of each simulation run for each grid cell and is plotted against the corresponding actual error (ϵ_u,v) in Fig. 1b. It seems that for this virtually steady case, there is a clear correlation between the diagnosed (instantaneous) χ values and the ϵ values that have accumulated over the simulation run time. Even though the correlation is not perfect, it provides a convenient and consistent framework for a grid-adaptation algorithm. As such, a second convergence test for this case is performed using a variable-resolution grid within the domain. The mesh is based on the aforementioned adaptive-grid approach. For these runs, we iteratively decrease the so-called refinement criterion (ζ_u,v) by factors of 2 between the runs and monitor the increasing fidelity of the numerically obtained solution for all runs. The refinement criterion presets a threshold value (ζ) for the estimated error χ that defines when a cell should be refined (χ>ζ) or alternatively, when it may be coarsened ($\mathit{\chi }<\mathrm{2}\mathit{\zeta }/\mathrm{3}$). Figure 2a presents the results and compares the local grid resolution used against ϵ_u,v for the various (colour-coded) runs. It appears that for all separate runs, the algorithm employed a variable resolution mesh and that this has resulted in a smaller range of the local error in the solution (ϵ), as compared to the equidistant-grid cases. The local error in the solution is also compared to the wavelet-based estimated error in the representation of the solution fields in Fig. 2b. Compared to the results from the equidistant-grid approach as presented in Fig. 1b, the spread of the χ and ϵ values is relatively small for the separate runs when the adaptive-grid approach is used. The most prominent reason for the finite spread is that the error (ϵ) was diagnosed after 1000 time steps. This facilitated errors in the solution that arise in the solution at a specific location (with a large χ value) to “diffuse” over time towards regions where the solution remains to be characterized by a small χ value (not shown). Also, since u and v are coupled (due to the background rotation), local errors that arise in the solution for u “pollute” the v-component solution and vice versa. Furthermore, a spread is expected because the tree-grid structure only allows the resolution to vary by factors of 2 (Popinet, 2011).

Finally, the global convergence characteristics and the speed performance of the two approaches are studied. The global error (η) in the numerically obtained solution is evaluated as

$$\begin{array}{}\text{(24)}& {\displaystyle }\mathit{\eta }=\underset{z_{\mathrm{bottom}}}{\overset{z_{\mathrm{top}}}{\int }}\left({\mathit{ϵ}}_u+{\mathit{ϵ}}_v\right)\mathrm{d}z.\end{array}$$

In order to facilitate a comparison between the methods, we diagnose the number of grid cells (N) used for the adaptive-grid run. Figure 3a shows that for both approaches the error scales inversely proportional to the number of grid cells used to the second power (i.e. second-order spatial accuracy in 1-D). The adaptive grid results are more accurate than the results from the fixed-grid approach when employing the same number of grid cells. Figure 3b shows that for both approaches the required effort (i.e. measured here in wall-clock time) scales linearly with the number of grid cells, except for the runs that require less than 1∕10 of a second to perform. The plots reveals that per grid cell there is computational overhead for the adaptive-grid approach. These results show that the numerical solver used is well behaved.

The following sections are devoted to testing the adaptive-grid approach in a more applied SCM scenario, where the turbulent transport closures are applied (see Sect. 2) and the set-up is unsteady. Here, the quality of the adaptive-grid solution has to be assessed by comparing against reference results from other SCMs, large-eddy simulations and the present model running in equidistant-grid mode.

Figure 4 Time-averaged profiles over the ninth hour of the run according to the GABLS1 inter-comparison scenario. For (a) the horizontal wind components and (b) the potential temperature. Results are obtained with the present adaptive-grid SCM (coloured lines), the LES models ensemble (i.e mean±σ) from Beare et al. (2006) (grey-shaded areas), and the present SCM, employing an equidistant and static grid with a 6.25 m resolution (dashed lines). For z>250 m, the profiles have remained as they were initialized.

Figure 5 Evolution of (a) the vertical spatial-resolution distribution and (b) the total number of grid cells for the GABLS1 inter-comparison case.

3.2 GABLS1

The first GABLS inter-comparison case focusses on the representation of a stable boundary layer. Its scenario was inspired by the LES study of an ABL over the Arctic Sea by Kosović and Curry (2000). The results from the participating SCMs are summarized and discussed in Cuxart et al. (2006); for the LES inter-comparison study, the reader is referred to the work of Beare et al. (2006). The case prescribes the initial profiles for wind and temperature, a constant forcing for momentum corresponding to a geostrophic wind vector, ${\mathit{U}}_{\mathrm{geo}}=\left\{\mathrm{8},\mathrm{0}\right\}$ m s⁻¹, and the Coriolis parameter $f=\mathrm{1.39}\times {\mathrm{10}}^{-\mathrm{4}}$ s⁻¹. Furthermore, a fixed surface-cooling rate of 0.25 K h⁻¹ is applied, and ${\mathit{\theta }}_{\mathrm{v},\mathrm{ref}}=\mathrm{263.5}$ K. The model is run with a maximum resolution of 6.25 m and a domain height of 400 m. The maximum resolution corresponds to six levels of tree-grid refinement, where each possible coarser level corresponds to a factor of 2 increase in grid size.

Figure 6 Comparison of the results obtained with the adaptive-grid SCM and the participating models in the work of Svensson et al. (2011) for the vertical profiles of (a) the virtual potential temperature and (b) the wind-speed magnitude, for 14:00 local time on 23 October. Panel (c) the evolution of the 10 m wind speed (U_10 m) on 23 October. For the model abbreviations used in the legend, see Svensson et al. (2011). The different shades of grey in plot (c) indicate observations from different measurement devices; see Svensson et al. (2011) for the details.

Figure 7 Vertical profiles of the wind-speed magnitude U obtained with the adaptive-grid (in colour) and the fixed equidistant-grid (dashed) SCMs. The 12 plotted profiles are obtained for 24 October with an hourly interval, starting from 01:00 local time. Note that the profiles are shifted in order to distinguish between the different times (with vanishing wind at the surface). The profiles of U are constant with height for z>1200 m.

Figure 8 Evolution of (a) the vertical spatial resolution and (b) the total number of grid cells, for the GABLS2 inter-comparison case. Two full diurnal cycles, corresponding to 23 and 24 October 1999 (ranging from the labels 1:00:00 to 3:00:00 on the x axis).

Due to the idealizations in the case set-up with respect to the reality of the field observations, the model results were not compared to the experimental data (Cuxart et al., 2006). However, for the SCMs, a reference was found in the high-fidelity LES results that tended to agree well between the various models. The LES results therefore serve as a benchmark for the results obtained with the present model. This facilitates a straightforward testing of the formulations and implementations of the physical closures used before we continue our analysis towards the full diurnal cycle. Inspired by the analysis of Cuxart et al. (2006) and their Fig. 3, we compare our SCM results with the 6.25 m resolution LES ensemble results. We focus on the profiles for the wind components and potential temperature averaged over the ninth hour of the simulation in Fig. 4. We observe that the present SCM is in good agreement with the LES results and is able to capture the vertical structure of the ABL, including the low-level jet. The differences are only minor compared to the variations in the results presented in the aforementioned GABLS1 SCM reference paper.

Note that in general, results are of course sensitive to the closure chosen to parameterize the turbulent transport, in our case given by Eqs. (6) and (11). In order to separate between the numerical effects of using grid adaptivity and the chosen physical closures, we define an additional reference case in which we run an equidistant-grid SCM. This model run employs a fixed 6.25 m resolution (i.e with 64 cells), but otherwise identical closures and numerical formulations. That is, we have switched off the grid adaptivity and maintain the maximum resolution throughout the domain. We can observe that results between both SCMs are in good agreement but that minor deviations are present. These discrepancies are on the order of magnitude of the refinement criteria and can be reduced by choosing more stringent values, which would result in using more grid cells. The evolution of the adaptive-grid structure is shown in Fig. 5a. We see that a relatively high resolution is employed near the surface, i.e. in the logarithmic layer. Remarkably, without any a priori knowledge, the grid is refined at a height of 150 m $<z<$ 200 m as the so-called low-level jet develops, whereas the grid has remained coarse above the boundary layer where the grid resolution was reduced to be as coarse as 100 m. From Fig. 5b we learn that the number of grid cells varied between 11 and 24 over the course of the simulation run.

3.3 GABLS2

The second GABLS model inter-comparison case was designed to study the model representation of the ABL over the course of two consecutive diurnal cycles. The case is set up after the observations that were collected on 23 and 24 October 1999 during the CASES-99 field experiment in Leon, Kansas, USA (Poulos et al., 2002). The case prescribes idealized forcings for two consecutive days that were characterized by a strong diurnal cycle pattern. During these days, the ABL was relatively dry, there were few clouds, and ${\mathit{\theta }}_{\mathrm{v},\mathrm{ref}}=\mathrm{283.15}$ K. The details of the case are described in the work of Svensson et al. (2011), which was dedicated to the evaluation of the SCM results for the GABLS2 inter-comparison. Compared to the original case prescriptions, we choose a slightly higher domain size of z_top=4096 m (compared to 4000 m), so that a maximum resolution of 8 m corresponds to nine levels of refinement.

In this section we place our model output in the context of the results presented in the work of Svensson et al. (2011), which, apart from the SCM results, also includes the results from the LES by Kumar et al. (2010). To obtain their data we have used the so-called “data digitizer” of Rohatgi (2018). Inspired by the analysis of Svensson et al. (2011) and their Figs. 10 and 11, we compare our results for the wind-speed magnitude ($U=\Vert \mathit{u}\Vert$) and virtual potential temperature profiles at 14:00 local time on 23 October in Fig. 6a and b, respectively. Here, we see that the results obtained with the present SCM fall within the range of the results as were found with the selected models that participated in the original inter-comparison. These models also employed a first-order-style turbulence closure and have a lowest model-level height of less than 5 m. The present modelled virtual potential temperature (θ_v) shows a slight negative vertical gradient in the well-mixed layer. This is a feature related to the use of the local K-diffusion description for the turbulent transport (see Sect. 2 and the work of Holtslag and Boville, 1993). Figure 6c presents a time series of the 10 m wind speed (U_10 m) during 23 October. Again the present model results compare well with the others. Next, in order to validate the grid-adaptivity independently from the closures used, we present the hourly evolution of the wind speed on 24 October against the results obtained with adaptivity switched off, using 512 equally spaced grid points in Fig. 7. A nearly identical evolution of the wind-speed profiles is observed and even the small-scale features in the inversion layer (i.e. z≈800 m) are present in the adaptive grid-model calculations. The corresponding evolution of the adaptive-grid structure is presented in Fig. 8, where the colours in the resolution plot appear to sketch a “Stullian” image, showing a prototypical diurnal evolution of the ABL (Stull, 1988). Apparently, the grid-adaptation algorithm has identified (!) the “surface layer” within the convective boundary layer (CBL), the stable boundary layer, the entrainment zone, and the inversion layer as the dynamic regions that require a high-resolution mesh. Conversely, the well-mixed layer within the CBL, the residual layer, and the free-troposphere are evaluated on a coarser mesh. The total number of grid cells varied between 24 and 44.

4 Discussion and conclusions

In this work we have presented a one-dimensional (1-D) single-column model (SCM) that employs a mesh whose resolution is varied adaptively based on the evolution of the numerically obtained solution. This is an attractive feature because it is a prerequisite to enable the computational effort required for the evaluation of numerical solution to scale with the complexity of the studied physical system. The adaptation algorithm based on limiting discretization errors appears to function very well for a wide variety of geophysical applications, e.g. 3-D atmospheric turbulence-resolving models (Van Hooft et al., 2018), tsunami and ocean-wave modelling (Beetham et al., 2016; Marivela-Colmenarejo, 2017; Popinet, 2011), hydrology (Kirstetter et al., 2016), two-phase microphysics (Howland et al., 2016), flow of granular media (Zhou et al., 2017), and shock-wave formation (Eggers et al., 2017). For these studies on highly dynamical systems, the adaptive-grid approach is chosen because it offers a more computationally efficient approach as compared to the use of static grids.

The present work does not include an in-depth assessment and discussion on the performance of the presented methods in relation to the computational speed. Even though this is an important motivation for the application of the adaptive-grid strategy to GCMs, the authors argue that an SCM is not suitable for speed-performance testing: the speed of single-column calculations is virtually never a critical issue. Only in 3-D mode, when SCMs are “stitched together” to enable the resolving of global circulations does the model's computational efficiency become an issue. Furthermore, the performance of an SCM that employs a few tens of cells is not a good indicator for the performance of a GCM that can employ billions of grid cells. For the latter, parallel computation overhead and the so-called memory bottleneck are important aspects. In contrast, for the SCM case, the complete instruction set and solution data can typically be loaded onto the cache memory of a single CPU's core. Nevertheless, for the readers' reference, the required run times for the different SCM set-ups presented herein are listed in Table 1, and Fig. 3b also presents quantitative results on this topic and shows that the adaptive-grid solver is well behaved.

Following the turbulence-resolving study of Van Hooft et al. (2018), the results presented herein are a proof of concept for future 3-D modelling, using RANS techniques. The authors of this work realize that the present SCM is a far cry from a complete global model and that more research and development is required before the method can be compared on a global-circulation scale. As shown by, e.g., Jablonowski (2004), a 3-D adaptive grid also allows a variable grid resolution in the horizontal directions. This further enables the computational resources to focus on the most challenging atmospheric processes where there is temporal and spatial variation in the horizontal-resolution requirement of the grid. Examples include the contrasting dynamics between relatively calm centres of high-pressure circulations and those characterizing stormy low-pressure cells. Also, in the case of a sea breeze event (Arrillaga et al., 2016), it would be beneficial to temporarily increase the horizontal resolution near the land–sea interface. As such, we encourage the use of this technique for those meteorologically challenging scenarios.