the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

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

### J. Antoon Hooft

### Stéphane Popinet

### Bas J. H. Wiel

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.

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.

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

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

with *Ri*_{b} the surface bulk Richardson number, which is defined
as

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},\phantom{\rule{0.125em}{0ex}}\mathrm{ref}}({\mathit{\theta}}_{v}-{\mathit{\theta}}_{\mathrm{v},\phantom{\rule{0.125em}{0ex}}\mathrm{ref}})$ . The virtual potential temperature is related to
the potential temperature (*θ*) and specific humidity (*q*) according
to

with ${R}_{\mathrm{v}}\phantom{\rule{0.125em}{0ex}}/\phantom{\rule{0.125em}{0ex}}{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

with *k*=0.4 the von Kármán constant, *z*_{1} 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

which conclude the description of the surface fluxes. The vertical flux
($\stackrel{\mathrm{\u203e}}{{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

where *K* is the so-called eddy diffusivity

*l* represents an effective mixing length,

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,

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

i.e. based on the gradient Richardson number,

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,

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*^{n} and *s*^{n+1}, this can be written

Rearranging the terms we write

to obtain a Poisson–Helmholtz equation for *s*^{n+1}, using the eddy
diffusivity calculated from the solution *s*^{n} (*K*^{n}). 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*,

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

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

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

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.

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.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
(−**∇**_{h}*P*) 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;

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}.

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*^{n}, *v*^{n}) within each grid cell is diagnosed and is
defined here as

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*^{n} is the value of the numerical solution within the
cell. Note that *a*^{n} 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

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.

## 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^{−1}, and the Coriolis parameter $f=\mathrm{1.39}\times {\mathrm{10}}^{-\mathrm{4}}$ s^{−1}. Furthermore, a fixed surface-cooling rate of
0.25 K h^{−1} is applied, and ${\mathit{\theta}}_{\mathrm{v},\phantom{\rule{0.125em}{0ex}}\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.

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},\phantom{\rule{0.125em}{0ex}}\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.

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.

Basilisk is a freely available (GPLv3), multi-purpose tool to solve partial differential equations and has its own website: http://www.basilisk.fr. The code contains solvers for Saint-Venant problems, the Navier–Stokes equations, electro hydrodynamics, and more; see http://basilisk.fr/src/README. A selection of examples can be viewed online: http://www.basilisk.fr/src/examples. The website also provides general information including installation instructions and a tutorial. Furthermore, for the work herein, interested readers can visit the model set-up pages, and links to their online locations are presented in Table 1. The data can easily be generated by running the scripts. Finally, a snapshot of the code used, as it was used in this the work, is made available via ZENODO, with doi link: https://doi.org/10.5281/zenodo.1203631.

All authors contributed to the content of this paper, and it should be viewed as a fruit of the many discussions we have had. Furthermore, SP wrote the Basilisk code; he also designed and implemented the grid-adaptation algorithm. The numerical experiments were set up and performed by JAvH. The writing was led by JAvH and organized in an iterative procedure with BJHvdW.

The authors declare that they have no conflict of interest.

The authors gratefully acknowledge the funding by the ERC Consolidator grant
(648666). The LES ensemble results used for the GABLS1 inter-comparison are
kindly made available by Bob Beare (online via
http://gabls.metoffice.com/), and we thank three anonymous reviewers
for their comments on the paper. Furthermore, we acknowledge that we are
indebted to those who contributed to the GNU/Linux project and/or the GNU
compiler collection (see http://www.gnu.org/).

Edited by: James R. Maddison

Reviewed by: three anonymous referees

Arrillaga, J. A., Yagüe, C., Sastre, M., and Román-Cascón, C.: A characterisation of sea-breeze events in the eastern Cantabrian coast (Spain) from observational data and WRF simulations, Atmos. Res., 181, 265–280, 2016. a

Baas, P., van de Wiel, B., van der Linden, S., and Bosveld, F.: From Near-Neutral to Strongly Stratified: Adequately Modelling the Clear-Sky Nocturnal Boundary Layer at Cabauw, Bound.-Lay. Meteorol., 166, 217–238, 2017. a

Beare, R. J., Macvean, M. K., Holtslag, A. A., Cuxart, J., Esau, I., Golaz, J.-C., Jimenez, M. A., Khairoutdinov, M., Kosovic, B., Lewellen, D., Lund, T. S., Lundquist, J. K., Mccabe, A., Moene, A. F., Noh, Y., Raasch, S., and Sullivan, P.: An intercomparison of large-eddy simulations of the stable boundary layer, Bound.-Lay. Meteorol., 118, 247–272, 2006. a, b

Beetham, E., Kench, P. S., O'Callaghan, J., and Popinet, S.: Wave transformation and shoreline water level on Funafuti Atoll, Tuvalu, J. Geophys. Res.-Oceans, 121, 311–326, 2016. a

Bosveld, F. C., Baas, P., Steeneveld, G.-J., Holtslag, A. A., Angevine, W. M., Bazile, E., de Bruijn, E. I., Deacu, D., Edwards, J. M., Ek, M., Larson, V. E., Pleim, J. E., Raschendorfer, M., and Svensson, G.: The third GABLS intercomparison case for evaluation studies of boundary-layer models. Part B: results and process understanding, Bound.-Lay. Meteorol., 152, 157–187, 2014. a

Boussinesq, J.: Théorie de l'écoulement tourbillonnant et tumultueux des liquides dans les lits rectilignes à grande section, vol. 1, Gauthier-Villars, 1897 (in French). a

Cuxart, J., Holtslag, A. A., Beare, R. J., Bazile, E., Beljaars, A., Cheng, A., Conangla, L., Ek, M., Freedman, F., Hamdi, R., Kerstein, A., Kitagawa, H., Lenderink, G., Lewellen, D., Mailhot, J., Mauritsen, T., Perov, V., Schayes, G., Steeneveld, G.-J., Svensson, G., Taylor, P., Weng, W., Wunsch, S., and Xu, K-M.: Single-column model intercomparison for a stably stratified atmospheric boundary layer, Bound.-Lay. Meteorol., 118, 273–303, 2006. a, b, c

Dunbar, T., Hanert, E., and Hogan, R.: A one-dimensional finite-element boundary-layer model with a vertical adaptive grid, Bound.-Lay. Meteorol., 128, 459–472, 2008. a

Eggers, J., Grava, T., Herrada, M., and Pitton, G.: Spatial structure of shock formation, J. Fluid Mech., 820, 208–231, 2017. a

Emauel, K. A.: Atmospheric convection, Oxford university press, 198 Madison avenue, New York, New York, USA, 1994. a

England, D. E. and McNider, R. T.: Stability functions based upon shear functions, Bound.-Lay. Meteorol., 74, 113–130, 1995. a

Grell, G. A., Peckham, S. E., Schmitz, R., McKeen, S. A., Frost, G., Skamarock, W. C., and Eder, B.: Fully coupled “online” chemistry within the WRF model, Atmos. Environ., 39, 6957–6975, 2005. a

Heus, T., van Heerwaarden, C. C., Jonker, H. J. J., Pier Siebesma, A., Axelsen, S., van den Dries, K., Geoffroy, O., Moene, A. F., Pino, D., de Roode, S. R., and Vilà-Guerau de Arellano, J.: Formulation of the Dutch Atmospheric Large-Eddy Simulation (DALES) and overview of its applications, Geosci. Model Dev., 3, 415–444, https://doi.org/10.5194/gmd-3-415-2010, 2010. a

Holtslag, A. and Boville, B.: Local versus nonlocal boundary-layer diffusion in a global climate model, J. Climate, 6, 1825–1842, 1993. a, b, c, d, e

Holtslag, A., Svensson, G., Baas, P., Basu, S., Beare, B., Beljaars, A., Bosveld, F., Cuxart, J., Lindvall, J., Steeneveld, G., Tjernström, M., and Van De Wiel, B. J. H.: Stable atmospheric boundary layers and diurnal cycles: challenges for weather and climate models, B. Am. Meteorol. Soc., 94, 1691–1706, 2013. a

Howland, C. J., Antkowiak, A., Castrejón-Pita, J. R., Howison, S. D., Oliver, J. M., Style, R. W., and Castrejón-Pita, A. A.: It's Harder to Splash on Soft Solids, Phys. Rev. Lett., 117, 184502, https://doi.org/10.1103/PhysRevLett.117.184502, 2016. a

Jablonowski, C.: Adaptive grids in weather and climate modeling, PhD thesis, University of Michigan, USA, 2004. a, b

Jablonowski, C., Oehmke, R. C., and Stout, Q. F.: Block-structured adaptive meshes and reduced grids for atmospheric general circulation models, Philos. T. Roy. Soc. A, 367, 4497–4522, 2009. a

Kirstetter, G., Hu, J., Delestre, O., Darboux, F., Lagrée, P.-Y., Popinet, S., Fullana, J.-M., and Josserand, C.: Modeling rain-driven overland flow: Empirical versus analytical friction terms in the shallow water approximation, J. Hydrol., 536, 1–9, 2016. a

Kosović, B. and Curry, J. A.: A large eddy simulation study of a quasi-steady, stably stratified atmospheric boundary layer, J. Atmos. Sci., 57, 1052–1068, 2000. a

Kumar, V., Svensson, G., Holtslag, A., Meneveau, C., and Parlange, M. B.: Impact of surface flux formulations and geostrophic forcing on large-eddy simulations of diurnal atmospheric boundary layer flow, J. Appl. Meteorol. Clim., 49, 1496–1516, 2010. a

Lenderink, G. and Holtslag, A. A.: An updated length-scale formulation for turbulent mixing in clear and cloudy boundary layers, Q. J. Roy. Meteor. Soc., 130, 3405–3427, 2004. a

Louis, J., Tiedtke, M., and Geleyn, J.: A short history of PBL parameterization at ECMWF, in: paper presented at the Workshop on Planetary Boundary Layer Parameterization, Eur. Cent. For Medium-Range Weather Forecasts, Reading, England, 1982. a

Marivela-Colmenarejo, R.: Numerical Perspective on Tsunami Hazards and Their Mitigation by Coastal Vegetation, PhD thesis, Virginia Tech, Blacksburg, Virginia, USA, 2017. a

Mellor, G. L. and Yamada, T.: Development of a turbulence closure model for geophysical fluid problems, Rev. Geophys., 20, 851–875, 1982. a

Neggers, R. A., Siebesma, A. P., and Heus, T.: Continuous single-column model evaluation at a permanent meteorological supersite, B. Am. Meteorol. Soc., 93, 1389–1400, 2012. a

Popinet, S.: Quadtree-adaptive tsunami modelling, Ocean Dynam., 61, 1261–1285, 2011. a, b, c, d, e, f

Popinet, S.: Time-implicit discretisation of reaction–diffusion equations, available at: http://www.basilisk.fr/src/diffusion.h (last access: 1 January 2018), 2017a. a

Popinet, S.: Multigrid Poisson-Helmholtz solvers, available at: http://www.basilisk.fr/src/poisson.h (last access: 1 January 2018), 2017b. a, b

Popinet, S.: Basilisk 18-02-16,

https://doi.org/10.5281/zenodo.1203631, 2018.

Poulos, G. S., Blumen, W., Fritts, D. C., Lundquist, J. K., Sun, J., Burns, S. P., Nappo, C., Banta, R., Newsom, R., Cuxart, J., Terradellas, E., Balsley, B., and Jensen, M.: CASES-99: A comprehensive investigation of the stable nocturnal boundary layer, B. Am. Meteorol. Soc., 83, 555–581, 2002. a

Reynolds, O.: On the dynamical theory of incompressible viscous fluids and the determination of the criterion, Philos. T. R. Soc. Lond., 56, 40–45, 1895. a

Rohatgi, A.: WebPlotDigitizer, available at: https://github.com/ankitrohatgi/WebPlotDigitizer, last access: 1 January 2018. a

Schaller, R. R.: Moore's law: past, present and future, IEEE Spectrum, 34, 52–59, 1997. a

Slingo, J.: The development and verification of a cloud prediction scheme for the ECMWF model, Q. J. Roy. Meteor. Soc., 113, 899–927, 1987. a

St-Cyr, A., Jablonowski, C., Dennis, J. M., Tufo, H. M., and Thomas, S. J.: A comparison of two shallow-water models with nonconforming adaptive grids, Mon. Weather Rev., 136, 1898–1922, 2008. a

Stull, R. B.: An introduction to boundary layer meteorology, vol. 1, Springer Science & Business Media, Springer, Dordrecht, The Netherlands, 670 pp., 1988. a

Svensson, G., Holtslag, A., Kumar, V., Mauritsen, T., Steeneveld, G., Angevine, W., Bazile, E., Beljaars, A., De Bruijn, E., Cheng, A., Conangla, L., Cuxart, J., Ek, M., Falk, M. J., Freedman, F., Kitagawa, H., Larson, V. E., Lock, A., Mailhot, J., Masson, V., Park, S., Pleim, J., Söderberg, S., Weng, W., and Zampieri, M.: Evaluation of the diurnal cycle in the atmospheric boundary layer over land as represented by a variety of single-column models: the second GABLS experiment, Bound.-Lay. Meteorol., 140, 177–206, 2011. a, b, c, d, e, f

Van Hooft, J. A., Popinet, S., Van Heerwaarden, C. C., Van der Linden, S. J. A., de Roode, S. R., and Van de Wiel, B. J. H.: Towards Adaptive Grids for Atmospheric Boundary Layer Simulations, Bound.-Lay. Meteorol., 167, 421–443, 2018. a, b, c, d, e, f

Wyant, M. C., Bretherton, C. S., Chlond, A., Griffin, B. M., Kitagawa, H., Lappen, C.-L., Larson, V. E., Lock, A., Park, S., De Roode, S. R., Uchida, J., Zhao, M., and Ackerman, A. S.: A single-column model intercomparison of a heavily drizzling stratocumulus-topped boundary layer, J. Geophys. Res., 112, D24204, https://doi.org/10.1029/2007JD008536, 2007. a

Zhou, Y., Lagrée, P.-Y., Popinet, S., Ruyer, P., and Aussillous, P.: Experiments on, and discrete and continuum simulations of, the discharge of granular media from silos with a lateral orifice, J. Fluid Mech., 829, 459–485, 2017. a