\documentclass[gmd, manuscript]{copernicus}



\usepackage{acro}
\DeclareAcronym{GCM}{
    short = GCM,
    long = General Circulation Model
}
\DeclareAcronym{AGCM}{
    short = AGCM,
    long = Atmospheric Global Climate Model
}
\DeclareAcronym{AOGCM}{
    short = AOGCM,
    long = Atmospheric Ocean Global Climate Model
}
\DeclareAcronym{IPCC}{
    short = IPCC,
    long = Intergovernmental Panel on Climate Change
}
\DeclareAcronym{ESM}{
    short = ESM,
    long = Earth System Models
}
\DeclareAcronym{GFDL}{
    short = GFDL,
    long = Geophysical Fluid Dynamic Laboratory
}
\DeclareAcronym{PBL}{
    short = PBL,
    long = Planetary Boundary Layer
}
\DeclareAcronym{KE}{
    short = KE,
    long = kinetic energy
}
\DeclareAcronym{ISA}{
    short = ISA,
    long = International Standard Atmosphere
}
\DeclareAcronym{3D}{
    short = 3D,
    long = three-dimensional
}
\DeclareAcronym{2D}{
    short = 2D,
    long = two-dimensional
}
\DeclareAcronym{DCMIP}{
    short = DCMIP,
    long = Dynamical Core Intercomparison Project
}
\DeclareAcronym{HS}{
    short = HS,
    long = Held \& Suarez
}


\begin{document}

\title{WAVETRISK adaptive dynamical core coupled to dry physics}
\Author[]{G}{Ching-Johnson}
\Author[2]{NK-R}{Kevlahan}
\Author[3]{T}{Dubos}
\Author[4]{F}{Hourdin}

\affil[1]{Department of Meteorology, University of Reading, Reading RG6 7BE, United Kingdom}
\affil[2]{Department of Mathematics and Statistics, McMaster University, Hamilton ON L8S 4K1, Canada}
\affil[3]{Laboratoire de M\'et\'erologie Dynamique, \'Ecole Polytechnique, 91128 Palaiseau, France}
\affil[4]{Laboratoire de M\'et\'eorologie Dynamique, Institut Pierre Simon Laplace, Paris, France}
\correspondence{Gabrielle Ching-Johnson (g.n.ching-johnson@pgr.reading.ac.uk)}

\runningtitle{Adaptive dry physics climate model}
\runningauthor{G Ching-Johnson, N Kevlahan, T Dubos, F Hourdin}

\received{}
\pubdiscuss{} 
\revised{}
\accepted{}
\published{}
\firstpage{1}

\def\fh#1{\textcolor{red}{#1}}

\maketitle

\begin{abstract}
General circulation models (GCMs) of the atmosphere couple a dynamical core for fluid dynamics that are primarily two-dimensional with a physics model for processes like radiation and convection that are primarily one-dimensional. In this paper we fully describe the dry physics model of~\cite{Frederic_these}, which includes a soil column, radiation, vertical diffusion and convective adjustment.  The climatology of the dry physics coupled with the {\sc wavetrisk-atmosphere} dynamical core is described.  We show that the dynamically adaptive version of {\sc wavetrisk-atmosphere} can be stably and accurately coupled to the dry physics without scale-aware parameterizations or modifications to the dynamics-based adaptivity criteria. The turbulence generated by {\sc wavetrisk}--dry physics GCM is characterized by longitude-latitude projections of vorticity, three-dimensional vorticity isosurfaces and vertical profiles of energy spectrum power laws, which vary between $-2.8$ in the boundary layer to $-5.2$ in the stratosphere for the rotational mode.  
\end{abstract}

\introduction  \label{intro}
The coupling of atmospheric and oceanic \acp{GCM} constitutes the backbone of global climate models (for which the acronym GCM is also used). Pioneered in the 1950's \citep{Phillips1956}, \acp{GCM} are specific numerical tools based on separating and coupling two fundamental components:
\begin{enumerate}
\item 
A ``dynamical core" that computes the ``large scale" fluid motions using approximate versions of the compressible or incompressible Navier--Stokes equations and transport of a thermodynamic prognostic variable, such as potential temperature or buoyancy.
\item  
A set of ``physics" models, called ``parameterizations'', that account for physics other than fluid dynamics (e.g., radiation, phase changes, subsurface thermal conduction and hydrology). The physics includes subgrid scale (SGS) models to approximate the effect of unresolved fluid motions on the resolved scales (e.g., turbulent diffusion).
\end{enumerate}
The parameterization of subgrid scale motions relies on Reynolds decomposition, which separates scales between a ``large scale", or average, resolved component one one side, and an unresolved turbulent component, or perturbation on the other. In the Reynolds framework, the large scale is defined formally as the ensemble average of random realizations of the full turbulent flow.
In this decomposition, the statistics of the SGS variables are assumed to be horizontally homogeneous and  ergodicity is invoked to assume that ensemble averages are equivalent to both horizontal space and time averages (either a running mean, defined in each point of the  continuous world, or a grid box average).
Similarly, the plane-parallel approximation is used for radiation. The radiation model assumes that all quantities are locally uniform horizontally, above a plane surface. For both SGS motions and radiation, the assumption of horizontal homogeneity results in physics models which are one-dimensional in the vertical.

\acp{GCM}, based on the separation of three-dimensional (fluid) dynamics and one-dimensional physics parameterizations, have proved to be very useful for understanding the climate machine, and in particular the link between convective and cloud processes and climate. From the beginning they have also been used successfully for prediction.
An early success, in the 1960s, was weather prediction over a time horizon of a few days. In the late 1970s climate models based on coupled atmospheric and oceanic \acp{GCM} successfully predicted the subsequent global warming. They played a great role in the raising concern about the risks of global warming. 
oday climate models are used systematically to model the climate at seasonal, decadal and centennial time scales. They are also used to evaluate the sensitivity of the climate to variations in input data, most importantly the level of greenhouse gas emissions. Although \acp{GCM} have reached a high degree of maturity and sophistication, many challenges remain.

Lastly, a major bottleneck of climate models is the enormous computational power and associated computational time required to produce useful results. Any increase in complexity, through the inclusion of more processes and higher resolution, requires more computational power. With only limited computational power, modellers are required to make numerous trade-offs between CPU time, physical and numerical accuracy, and the number of physical processes represented.

Choosing horizontal grid resolution $\Delta x$ is the fundamental trade-off when designing and running a climate model.  More resolution provides higher accuracy, but greatly increases the computational cost.  Since the time step $\Delta t\propto 1/\Delta x$, doubling the resolution from $\Delta x$ to $\Delta x/2$ means a simulation year takes 8 times longer to compute on a given CPU (assuming memory is not a limitation), and in any case 8 times more operations to be done !
This is a main motivation for adaptative grids, in which the resolution is refined locally depending on the structure of the flow. 
\fh{I added the sentence before because I find it difficult to introduce the idea of adaptivity in a parenthesis for the vertical grid without introducing first for the horizontal grid}.
Note that only the horizontal resolution is adapted usually even if Lagrangian vertical coordinates with adaptive remapping and dormant layers could in principle provide vertical resolution adaptivity. Often, high resolution is not required everywhere and so uniformly increasing resolution wastes computational resources: some regions are over-resolved at some times and others remain under-resolved. Dynamically adaptive climate simulations have been proposed as a way of using computational resources optimally and ensuring uniform accuracy~\citep[e.g.,][]{dynam_adapt_Jablonowski,Kevlahan/Dubos:2019}. A dynamically adaptive GCM locally refines and coarsens the grid as a function of time and location based on dynamical, physical or numerical accuracy criteria~\citep{dynam_adapt_Jablonowski}.  Beyond capturing the evolution of the small scale dynamics,  adaptivity aids in capturing the effects of the interactions between small-scale and large scale processes~\citep{dynam_adapt_Jablonowski}. 

In theory, adaptive grids could improve the realism and accuracy of model simulations, but there are known complications and challenges that come with adaptivity. Most notably, as mentioned above, SGS physics parameterizations are often thought for a particular horizontal resolution and may not perform properly when the grid is refined and coarsened during the simulation~\citep{Collins13_grid_adapt}. This is expected to be the case in particular for moist deep convection, when the grid resolution reaches a few tens of kilometres, or for dry or cloudy boundary layer convection for resolutions of 1 or 2~km \citep{Honnert2016}. The question of whether physics processes should be ``scale-aware" in this ``grey zones" of convection is one challenge that many researchers are currently exploring~\citep[e.g.,][]{ scale_aware2,scale_aware1}. 


Scale-aware physics parameterizations adjust themselves based on the current local resolution of the grid. At its most basic, a scale-aware physics parameterization adjusts its internal model appropriately in response to changes in local grid resolution. A simple example is resolution dependent grid-scale viscosity, where the non-dimensional viscosity $\nu\Delta t/\Delta x^2$ is kept constant when the grid resolution $\Delta x$ changes. In the most extreme case where the resolution is refined enough to allow for the process to be fully resolved, the SGS parameterization is deactivated. Specifically, in the case of dynamic adaptivity, a refinement criteria needs to be met before local refinement (or coarsening) can occur. A central challenge of adaptivity is whether these criteria needs to be based on both the physics and the dynamics, or where adapting on the dynamics alone is sufficient. The answer likely depends on the type of physics being modeled. In current models the criteria are dynamics-based. For example, based on gradients or vorticity \citep[e.g.,][]{FlowCriteria}, or based on the numerical approximation~\citep[e.g.,][]{ErrorCriteria, Adaptive_Mesh_criteria,Kevlahan/Dubos:2019}. In order to make progress addressing these questions, researchers need well-understood physics models that can work in both the adaptive and non-adaptive models.

One strategy to tackle a particular bottleneck and test new approaches such as adaptative grids, is indeed to simplify other components or aspects of the model in order to avoid the complexity of running and analyzing results of a full \acp{GCM}.
For example, models may be simplified by reducing the number or complexity of components included. \acp{AGCM} simulates the general circulation of the atmosphere, while an \ac{AOGCM} couples the ocean and atmosphere sub-models to simulate the general circulation.  The atmosphere and ocean \acp{GCM} can also be more or less complex, and include more or fewer physical processes. When addressing botelnecks that concern the atmospheric dynamics, the ``physics" can be simplified further. The most basic atmosphere \acp{GCM} are hydrostatic and neglect topography, moist physics and ocean--atmosphere coupling. The land surface model is reduced to a basic radiation boundary condition and turbulent boundary layer model. Such highly simplified models are useful to better understand physics--dynamics coupling and the effect of resolution.

A classical way new dynamical cores or approaches is to replace the atmospheric physics by a relaxation toward a prescribed latitude-altitude temperature field, aiming at replacing in particular both radiation and SBS transport, with drag in the first layer~\citep{HeldSuarez}.
\fh{Do you really need the sentences below ? What for ? The wording "aquaplanet" is often used to describe an atmospheric model forced by SSTs that depend on latitude only. What you describe here as an aquaplanet is rather known at atmospheric GCM with slab ocean. How are those citation related to the present work ? Do you present series of papers that were devoted to some aspects of the dynamical core but using various configurations of the physics ? If it is the idea, it should be saied more clearly I feel. Also, we should then say why we make a particular choice compared to those other studies}
These models can be classified as dry physics ~\citep[e.g.,][]{Schneider2004, SchneiderWalker2006, MbengueSchneider2013,MbengueWoolings2019,ChemistryDry}, moist physics ~\citep[e.g.,][]{moist_test_case} and aquaplanet physics schemes ~\citep[e.g.,][]{moist_dynam_cores,Frierson2006,OGorman2008,Merlis2013a,Merlis2013b,JuckerGerber,Clark2018,Clark2020}, all with a range of complexities. The dry and moist physics models are physics models that must be coupled to a dynamical core, while the aquaplanet models are idealized atmospheric climate models that include a slab of ocean serving as the planet's surface. 

In this paper, we make use of a simple and versatile dry physics model developed by \citet{Frederic_these} to explore various atmospheric general circulation regimes for various planetary conditions. This model includes sub-models for radiation, convection, turbulence, the atmospheric boundary layer and the soil. This simple physics was used to investigate, changing only the value of a few parameters, some aspects of the terrestrial atmospheric circulation (Haddley circulation and it sensitivity to planetary roration rate and radiative properties), as well as the ciruclation of the Mars, Venus and Titan. It was used to investigate in particular atmospheric ``superrotation'' \citep{Hour:92b}, a phenomenon dominates the circulation of Venus and Titan. 

We chose this particular dry physics model to explore physics--dynamics coupling in a fairly realistic, yet simple, dry physics model with all assumptions and approximations fully described. The lack of moisture keeps the physics simple, but is still complex enough for us to investigate basic physics-dynamics coupling questions.  By using a dry physics, without phase change, it pushes the scale at which parameterizations of SGS motions start to be questionable down to a few kilometers. A strength of this physics model is that it can be easily coupled with different dynamical cores and adapted to model different planets, including gas giants.  The final reason for focusing on this particular model is to fully document, implement and test a package that can be coupled to both adaptive and non-adaptive dynamical cores. This SGS physics parameterization package  allows investigation of dynamic adaptivity, scale-awareness and the scale-dependence of basic physics processes, and could help reduce the need of climate model tuning with increases in resolution. 

With these goals in mind, this paper presents the simple dry physics model in detail and evaluates its climatology in the non-adaptive and adaptive cases by coupling it to an existing dynamical core. This initial characterization of the simple physics model provides guidance on physics-dynamics coupling for adaptive GCMs, and introduces the simple physics package as a way of comparing different dynamical cores by coupling them to the same physics package.  Specifically, this paper addresses the following key questions:
%
\begin{enumerate}
\item 
What is the climatology of the {\sc wavetrisk}--simple physics climate model? How does it differ from the \cite{HeldSuarez} climate model?

\item 
Is the dynamically adaptive \ac{GCM} numerically stable and and able to capture the emergence of small scale structure? Is adaptivity efficient for seasonal climate simulations?

\item 
Is scale-aware dry physics and/or physics dependent adaptivity required for stable and accurate dynamically adaptive climate simulations? 
\item 
Does the adaptive simulation capture the full range of active turbulent length scale in the atmosphere, as well as the correct power law scaling?

\item 
Can the {\sc wavetrisk}--simple physics climate model be restarted successfully from a checkpoint with significantly increased resolution?
\end{enumerate}

The paper is organized as follows. Section~\ref{SimplePhysics} describes in detail the simplified dry physics model, explaining its components and discussing its assumptions, limiting features and potential improvements. Section~\ref{Configuration} describes the coupling of the dry physics model with the adaptive dynamical core {\sc wavetrisk}~\citep{wavetrisk,Kevlahan/Lemarie:2022} in both the adaptive and non-adaptive cases. Section~\ref{Results} presents and interprets the results of the physics-dynamics coupling in the non-adaptive and adaptive cases, exploring the climatology and the effects of seasons, and Section~\ref{Conclusions} summarizes our conclusions and proposes future research directions.

\section{Simple Dry Physics Model}\label{SimplePhysics}
This section introduces the simple dry physics model.  We begin by summarizing the model's assumptions and then describe in detail each component of the model (radiation, turbulent diffusion, convection, planetary boundary layer, soil column), with a particular focus on boundary conditions and time integration. Lastly we describe the limitations of the package and its potential improvements. 
As with most physics models, the simple dry physics is implemented as a time integration split step and works on a set of vertical columns (of varying cross-section in the adaptive case). The fact that the dynamics works on individual horizontal layers and the physics works on individual vertical columns is a fundamental approximation of GCMs, and is well-suited to horizontal grid adaptivity.

\subsection{Physics model assumptions}

To properly understand a physics model it is essential to explicitly state its assumptions and approximations. This not only allows for full transparency, but also avoids inconsistent assumptions when the physics and dynamics are coupled~\citep{Lauritzen/etal:2022}. 

The basic assumptions of this physics package are: dry thermodynamics, compressible flow (i.e.\ atmosphere, not ocean),  hydrostatic pressure, and a floating (Lagrangian) vertical hybrid coordinate.

{\bf No moisture and compressibility approximations.}
A fundamental approximation of this package is the assumption that moisture (e.g., water in its three phases) is not taken into account. In particular, condensible substances and their effects on the atmospheric dynamics are absent. Such models are called dry physics models. Secondly, the model assumes compressible flow, which means it is suitable for the atmosphere, not the ocean. A fluid is compressible if changes in pressure causes significant changes in density. When modeling the atmosphere, the density of dry air changes significantly with changes in atmospheric pressure. Therefore, these two basic assumptions of no moisture and compressibility need to be ensured when coupling with a dynamical core. 

{\bf Hydrostatic approximation.} Another fundamental approximation in this simple physics model is the hydrostatic vertical pressure. The coupled dynamical model must therefore also assume hydrostatic balance.  The hydrostatic approximation constrains the vertical coordinates in the model and simplifies the vertical equation of motion. An important implication of the hydrostatic approximation is that dry convection must be parameterized and cannot be resolved (even if the local grid resolution permits it).

{\bf Floating vertical coordinate.}
There are many possible spatial coordinate systems that can be used when modeling a planet. For example, the Cartesian coordinate system, with coordinates $(x, y , z)$, where $z$ is the vertical height above mean sea level. In the case of simplified \acp{AGCM}, where the dynamics and the physics are coupled, the dynamics works on two-dimensional layers and the physics works on one-dimensional vertical columns. While $z$ might be one's first choice for a vertical coordinate, other possible vertical coordinates include pressure, mass and sigma (a normalized pressure coordinate).  

Vertical fluxes of heat causes air parcels to expand and compress, and in turn change the density $\rho$ (and associated buoyancy), which is represented in the mass continuity equation by
%
\begin{equation}
    \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial z}(w\rho) = 0,
\end{equation}
%
\noindent where $w$ is the vertical velocity.

When altitude $z$ is the vertical coordinate there are vertical fluxes across the horizontal layers and the dynamics must include vertical velocity as a prognostic variable (in addition to horizontal velocities in each layer). The vertical velocity must be diagnosed from the density, which can be a tedious process. An easier option, implemented in this physics package, is to use a Lagrangian vertical coordinate, sometimes known as a floating vertical coordinate, so there is no mass flux across horizontal layers. In the case of a Lagrangian vertical coordinate $\zeta$, the $z$ coordinate has the form $z = z(\zeta,t)$ and there is no explicit vertical velocity (and no mass flux across vertical layer interfaces) since the horizontal layers move vertically with the local vertical velocity. Therefore, the incremental mass $dM$ is
%
\begin{equation}
    dM = \rho dz = m\, d\zeta = \rho\frac{\partial z}{\partial \zeta}\partial \zeta
\end{equation}
%
and $dM/dt = 0$, where $m = \partial M/\partial \zeta  = \rho\partial z/\partial \zeta$.
%
Note that since the thickness of the horizontal layers changes over time, the vertical coordinates must be remapped periodically onto the initial vertical coordinates, or optimized target vertical coordinates.

\subsection{Model parameterizations}
%
The physics package includes parameterizations for radiation, vertical turbulent diffusion of velocity and potential temperature, a planetary boundary layer surface scheme, dry convective adjustment, soil column (for realistic surface flux boundary conditions), and diurnal radiation forcing. All models are one-dimensional (vertical) in space, and therefore are applied independently to each one-dimensional vertical column on the sphere as an implicit split step after the dynamics time step (i.e., the physics package takes as input a set of vectors). Note that none of the physics sub-models explicitly require the horizontal cross-section of the column (i.e, its horizontal resolution), although horizontal grid resolution could implicitly modify the physics models and physics--dynamics coupling. This is one of the questions we are interested in. 

The time scheme for the physics step uses Strang splitting:
\begin{enumerate}
\item 
Explicit forward Euler for soil, surface flux, radiation, and dry convective adjustment.
\item 
Implicit backwards Euler for vertical turbulent diffusion of velocity and potential temperature.        
\end{enumerate}
In the following, we give full details of each sub-model, including their default parameters. More details on the physics time scheme is given in Section~\ref{sec:physics_time}.

\subsubsection{Radiation}
%
Radiation is a process that has a major impact on the atmosphere. In general, since radiation is not a dynamical process and due to its complexity, it's parameterization involves a hierarchy of simplifying assumptions (e.g., local thermodynamic equilibrium, only vertical radiative transfer, simplified or neglected scattering, discretized spectral bands, \ldots). 

The radiation model included in the physics package is based on the two-stream approximation. The two-stream approximation equations make the plane-parallel assumption, in which radiation is a function of height and the atmosphere is horizontally homogeneous and stratified in flat, horizontal layers (at least on the horizontal scale of the column). The resulting downward and upward radiation fluxes at a specific frequency $\nu$ and height $z$ are 
%
\begin{align}
    F_{+}(\nu, z) &= F_+(\nu, 0) \tau(\nu,0, z) +  \int_{0}^{z'} \pi B(\nu,T(z')) K(\nu,z')\tau(\nu,z', z) \,dz' \nonumber\\ 
    &= F_+(\nu, 0) \tau(\nu,0, z) +  \int_{0}^{z'} \pi B(\nu,T(z')) \frac{\partial\tau(\nu,z', z)}{\partial z'} \,dz' \\
    F_{-}(\nu, z) &= F_-(\nu, \infty) \tau(\nu,  z, \infty) +  \int_{z'}^{\infty} \pi B(\nu,T(z')) K(\nu,z')\tau(\nu,z, z') \,dz' \nonumber\\
    &= F_-(\nu, \infty) \tau(\nu,z, \infty) -  \int_{z'}^{\infty} \pi B(\nu,T(z')) \frac{\partial\tau(\nu,z, z')}{\partial z'} \,dz',
\end{align}\label{two-stream}
%
where the $+$ subscript indicates an upward flux and the $-$ subscript indicates a downward flux. $B(\nu,T(z')$ is Planck's function: 
%
\begin{equation}
    B(\nu,T(z')) = \frac{2h\nu^3}{c^2(e^{h\nu /kT(z')}-1)},
\end{equation}
%
where $h$ is Planck's constant. Planck's function describes the radiation emitted by a black body in thermal equilibrium at a given temperature $T(z')
$.

$\tau_{\nu}(z, z')$ is the transmission function, which represents the proportion of the incident radiation, at frequency $\nu$, which is emitted from layer $z$ to layer $z'$. Therefore, it is the proportion of the flux that is not absorbed and in this model it is crudely represented by
%
\begin{equation}
    \tau(z, z') = e^{-|\psi(z) - \psi(z')|^\alpha},
\end{equation}
%
where $\psi$ is a monotonic function and $\alpha$ represents the absorption weight. This is known as a band averaged transmission function, as it approximates the transmission function over a given frequency band. The particular transmission function differs depending on the radiation band, thus the monotonic function in the physics package depends on the frequency band. In general, transmission functions are multiplicative. However, due to band averaging of the transmission function this property is lost.  

$K(\nu,z')$ represents the absorption, at layer $z'$, of radiation at frequency $\nu$, and is dependent on the infrared absorbers (i.e.\ gas) in the atmosphere,
%
\begin{equation}
    K(\nu,z') = \sum_{i}{}  \rho_i \kappa_i(p, T, \nu),
\end{equation}
%
where $\kappa(p, T, \nu)$ is the absorption coefficient of gas $i$ and $\rho$ is the density of the gas $i$.  The simple physics model assumes that the absorption coefficients depend only on height.

The simple physics model assumes constant absorption coefficients for the longwave (visible) and shortwave (visible) bands. The longwave (infrared) absorption coefficient is $\kappa_L = 0.08$\,\unit{m^{-1}} and the shortwave (visible) absorption  coefficient is $\kappa_S =0.99$\,\unit{m^{-1}}  . 

In general, the transmission function $\tau_{\nu}$ is a function of $K(\nu,z')$. Due to the band averaging of the transmission function, $K(\nu,z')$ is not explicitly defined. However the band averaged transmission function does describe strong or weak absorption of a frequency band through $\alpha$. A strong absorption produces a steep decay of transmission and strong absorption of radiation over shorter paths~\citep{pierrehumbert_2010}. This means that for longer paths the absorption in the first part of the path is stronger, while a  absorption is weaker further along in the path. Therefore, strong absorption is usually used for frequencies that are absorbed first in a path.

To solve for the net upward and downward irradiance, the fluxes need to be integrated over all frequencies. The physics packages uses the dual band two-stream approximations by crudely integrating the fluxes over just two spectral bands, longwave (L) and shortwave (S). Therefore, the net flux of radiation is
%
\begin{equation}
    F_{rad} = \int_{0}^{\infty} F_-(\nu,z) - F_+(\nu,z) \,d\nu = F_L(z) + F_S(z),
\end{equation}
%
where $F_L$ and $F_S$ are respectively the net longwave and shortwave radiation fluxes. 

\subsubsection{Shortwave radiation}
%
The shortwave net flux is the net result of the upward and downward fluxes evaluated for shortwave radiation,
%
\begin{equation}
    F_S(z) =  F_{S+}(z) - F_{S-}(z).
\end{equation}
%

The simple physics case neglects atmospheric emission of shortwave radiation, leaving the single initial terms as the upward and downward flux,
%
\begin{align}
    F_{S-}(z) &= \tau_S(z,\infty) \mu S \label{downShort}\\ 
    F_{S+}(z) &= \tau_S(0,z)F_{S+}(0) \label{upShort}
\end{align}
%
where $S$ is the downward flux at the top of the atmosphere, which in the case of shortwave radiation is the solar flux and $\mu$ is the cosine of the zenith angle. The solar flux represents the amount of flux from the  solar star impinging on the outer atmosphere. The upward flux at height $z$ (Equation \ref{upShort}) represents the proportion of the upward flux at the ground transmitted through the atmosphere, while the downward flux represents the proportion of the incident solar flux transmitted through the layer at $z$. 

The boundary upward flux of radiation, $F_{S+}(0)$, is the portion of downward shortwave radiation that is reflected by the ground,
%
\begin{equation}
    F_{S+}(0) = -(1-\alpha)F_{S-}(0),
\end{equation}
%
where $\alpha$ is the albedo, which depends on the reflective properties of the ground ($\alpha\approx 1$ for fresh snow, $\alpha\ll 1$ for forests). The albedo $\alpha=0.112$ for our simulations, independent of position, although the physics package allows location-dependent albedo (e.g., different for land and ocean regions). Therefore there is no shortwave radiation emitted by the ground itself, however the ground does emit longwave (heat) radiation and is taken into consideration in the longwave fluxes described below.

The transmission factor $\tau_S$ for the shortwave spectral band depends on height according to 
%
\begin{align}
    \tau_S(z,\infty) &= e^{-c_S \frac{\psi(z)}{\mu}}, \nonumber\\
    \tau_S(0,z) &= e^{-c_S\frac{\psi(0)-\psi(z)}{\mu_0}}, \nonumber\\
    \psi(z) &= \frac{p(z)}{p_{\mathrm rad}} \label{eq:trans_factor},
\end{align}
%
where $p_{\mathrm rad}=1000$\,\unit{hPa} is a reference pressure height and $c_S = 0.005$ is the attenuation of shortwave radiation at pressure height $p_{\mathrm rad}$. $\mu$ is the cosine of the zenith angle and $\mu_0$ is an average of $\mu$,  $\mu_0=3/5$. 

\subsubsection{Longwave radiation}
%
The longwave radiative fluxes in the simple physics package are represented using the gray gas approximation of the two-stream equations. As in the case of the shortwave net flux, the longwave net flux is determined by 
%
\begin{equation}
    F_L(z) =  F_{L+}(z) - F_{L-}(z), 
\end{equation}
%
where $F_{L+}(z)$ and $F_{L-}(z)$ represent the upward and downward longwave fluxes respectively. 

A further assumption in the gray gas model is the neglect of longwave solar emission meaning $F_{L-}(\infty) = 0$. Also, since the longwave radiative fluxes are the result of the two-stream approximations (\ref{two-stream}) integrated over infrared frequencies, the integration over infrared frequencies of the second term is well approximated by the Stefan--Boltzmann law~\citep{pierrehumbert_2010}. More specifically,
%
\begin{equation}
    \int \pi B \frac{\partial}{\partial z'}\tau(\nu,z,z') \,d\nu \approx \sigma T(z)^4 \frac{\partial}{\partial z'}\tau_{L}(z,z'),
\end{equation}
%
where $\sigma=5.67\times 10^{-8}$~W m$^{-2}$~K$^{-4}$ is the Stefan--Boltzmann constant. Therefore, the downward and upward flux of longwave radiation can be approximated as
%
\begin{align}
    F_{L-}(z) &= \int F_{-} \,d\nu = -\int \sigma T(z')^4\frac{\partial \tau_L(z,z')}{\partial z'}dz'\\
    F_{L+}(z) &= F_{L+}(z=0)\tau_L(0,z) + \int \sigma T(z')^4\frac{\partial \tau_L(z',z)}{\partial z'}dz'.
\end{align}
%
The transmission factor $\tau_L$ for the shortwave spectral band depends on height according to (\ref{eq:trans_factor}), but with strong absorption coefficient $C_L= 2.53$. 

The upward surface boundary radiation flux is comprised of the heat flux emitted from the ground itself and the fraction of the downward atmospheric longwave radiation that is reflected by the ground (i.e., not absorbed by the ground). Therefore,
%
\begin{equation}\label{boundaryuprad}
    F_{L+}(0) = \epsilon \sigma T_s^4 + (1-\epsilon)F_{L-}(0),
\end{equation}
%
where $\epsilon = 1 $ and $T_s$ are the soil emissivity and the surface temperature. The emissivity and the albedo, mentioned in the shortwave section are both characteristics of the ground, but they describe the effects for different wavelengths. The albedo $\alpha$ is the fraction of shortwave radiation reflected by a surface, while the emissivity is the effectiveness of a surface, the ground in this case, to emit longwave (heat) radiation. Both parameters are an input to the initialization of the physics package. Rearranging  (\ref{boundaryuprad}) gives the net longwave flux at the boundary,
%
\begin{equation}
    F_{L}(0) = F_{L+}(0) - F_{L-}(0) = \epsilon (\sigma T_s^4 - F_{L-}(0)).
\end{equation}

The transmission factors for the longwave spectral band are represented as
%
\begin{align}
    \tau_L(z,z') &= e^{-c_{LW} \sqrt{|\psi(z')-\psi(z)|}}\\
    \psi(z) &=\left(\frac{p(z)}{p_{\mathrm rad}}\right)^2,
\end{align}
%
where $c_{LW} = 2.53$ is the attenuation of the longwave radiation at pressure height $p_{\mathrm rad}=1000$\,\unit{hPa}. (A weak absorption option is also available for longwave radiation with $\psi(z) = p(z)/p_{\mathrm rad}^2$ .)

Upon discretization, the longwave equations are computed as,
%
\begin{align}
    F_{L-}(l) &=  \sum_{k>l}{} (\tau_L(l,k_-) - \tau_L(l,k_+))\sigma T(k)^4\\
    F_{L+}(l) &= F_{L+}(z=0)\tau_L(0,z) + \sum_{k<l}{} (\tau_L(k_+,l) - \tau_L(k_-,l))\sigma T(k)^4,
\end{align}
%
where $l$ is the layer interface and $k$ is the center of the layer. Therefore $k_+$ and $k_-$ represent the interfaces above and the below layer $k$ respectively.

The warming or cooling of a layer is affected by the net change in radiation with height, and is derived from the conservation of energy equation. This gives an expression for the temperature tendency,
%
\begin{equation}
    \frac{\partial T}{\partial t} = -\frac{1}{\rho c_p}\frac{\partial F_{rad}}{\partial z},
\end{equation}
%
where $c_p=1004$\,\unit{J\,kg^{-1}\,K^{-1}} is the specific heat capacity of air~\citep{IntroMicrometeorology}. 

Finally, the change in net radiative surface boundary flux with respect to the surface temperature is
%
\begin{equation}
    \frac{\partial F_{rad}(l)}{\partial T_s} = 4 \epsilon \sigma T_s^3\tau_L(0,l).
\end{equation}
%
This relation is used for the implicit coupling at the boundary.

\subsubsection{Astronomical parameters and diurnal cycle}
% 
The radiative transfer model is forced by solar radiation, which is time-dependent and set by the planet's astronomical parameters (e.g., obliquity, which determines seasons) and the diurnal (day/night) cycle.  In this paper we set these parameters to their Earth values, which are listed in Table~\ref{tab:astronomical}. The physics package includes a flag that activates diurnal radiation forcing.

\begin{table}
\begin{tabular}{|c|c|} 
\hline
{\bf Parameter} & {\bf Value} \\ \hline
Solar constant & $1370$\,\unit{W\,m^{-2}} \\
Perihelion distance & $150\times 10^6$\,\unit{km} \\
Ellipticity & 0 \\
Orbital period & 1 year \\
Planet radius & 6371\,\unit{km} \\
Gravitational acceleration $g$ & 9.8\,\unit{m\,s^{-2}} \\
Reference temperature $T_0$ & 285\,\unit{K} \\
Obliquity & $23.5^\circ$ \\
Rotation $\Omega$ & $7.292\times 10^{-5}$\,\unit{rad\,s^{-1}} \\
Day length & 24\,\unit{h} \\ \hline
\end{tabular}
\medskip
\caption{\label{tab:astronomical} Astronomical parameters used for the simple physics package (based on Earth values).}
\end{table}

\subsubsection{Vertical turbulent diffusion}
%
Turbulence has a major impact on the dynamics of atmospheric flow. In general, turbulence spans  a large range of spatial and time scales, from sub-millimetre to kilometres. More importantly for us, the computational cost of computing an integral time scale of a homogeneous three-dimensional turbulent flow scales like $Re^3$, where $Re$ is the non-dimensional Reynolds number characterizing the turbulence. Since $Re=O(10^8)$ for the atmosphere, it is clear that turbulence cannot be fully represented at the resolved scales and at least some turbulent length scales must be parameterized. This parameterization represents the effects of turbulence at the unresolved small scales on the resolved large scales. 

In GCMs turbulence, like other dynamics and physics, is decomposed into its horizontal and vertical components.  In our simulations, and climate simulations in general, the inertial (power law energy spectrum) range of horizontal turbulence is partially resolved.  Since the minimum resolution is in the inertial range the unresolved horizontal turbulence needs to be modeled.  In our case, we use second-order hyperdiffusion (i.e. $\nu_2 \Delta^2 u$) to model the effect of unresolved turbulence on the resolved turbulence as diffusion.  The kinematic viscosity $\nu_2$ is chosen to be the minimum required to remove grid-scale noise based on the non-dimensional viscosity. In the adaptive case, diffusion can either be based on the maximum grid resolution, or it can be scale-aware (based on the local grid resolution).  It is also possible to use more sophisticated horizontal turbulence models, such as large eddy simulation (LES).

In contrast, GCMs never resolve any part of the the vertical turbulence in each column: the effect of vertical turbulence is modeled as an enhanced vertical diffusion of velocity and potential temperature.  The goal of these models is therefore to approximate the vertical turbulent viscosity for velocity and temperature (i.e., eddy viscosity, turbulent diffusivity) for each column based on inputs from the dynamics step.

The simplest and most common model for vertical turbulent viscosity is the first-order one equation mixing length closure approximation. Ocean models such as NEMO use the second-order closure $k$-$\epsilon$ model, which solves a coupled dynamical system for vertical turbulent kinetic energy (TKE) $k$ and turbulent dissipation $\epsilon$ in order to compute the the effective turbulent viscosity coefficients.

Recall that the prognostic variables are assumed to be quantities averaged over large scales $L$ and long times $T$ (i.e.\ they are the ``large scale'' or coarse-grained resolved by the dynamics model).   The effect of the unresolved vertical turbulence on the prognostic variables in each layer is modeled as an enhanced vertical diffusion in their tendencies. The vertical turbulent diffusion model therefore converts a portion of the horizontal TKE into heat at a rate $\epsilon$ that depends on the diffusion coefficient $K_z^u$ and the vertical shear $\partial_z \textbf{u}$. To ensure conservation of energy, $\epsilon$ must be included as a source term in the heat equation~(\ref{turb_heat}).  Vertical turbulent diffusion also transfers velocity and potential temperature between horizontal layers. 

Since turbulence is assumed to act as an enhanced diffusion, the vertical flux of a quantity $q$ is given by 
%
\begin{equation} \label{turbflux}
    F_q = - \rho K_z^q \frac{\partial q}{\partial \zeta},
\end{equation}
%
where $\rho$ is the density, $K_z^q$ is a variable turbulent viscosity (to be modeled), and $\partial q/\partial \zeta$ is the vertical gradient of $q$. 

In general, the turbulent viscosity $K_z^q$ depends on the variable $q$. However, the physics package assumes that this coefficient is the same for both velocity and potential temperature so eddy diffusivity and  eddy viscosity are assumed to have the same values: $K_z^\theta = K_z^u$~\citep{Frederic_these}.

The eddy viscosity is computed based on a mixing length model, where $K_z^u$ is proportional to the vertical shear and the square of a ``mixing length'' $l$. The mixing length is the mean distance over which the quantity $q$ becomes fully mixed by the turbulence.  In addition, the eddy viscosity is proportional to a nonlinear term that measures the net production of turbulence by buoyancy and velocity shear instabilities.
%
\begin{equation} 
    K_z^u = l(z)^2 \left|\left|\frac{\partial \textbf{u}}{\partial z}\right|\right| \sqrt{1-\frac{Ri}{Ri_c}}.
    \label{eq:eddy_viscosity}
\end{equation}
%
Computing the eddy viscosity $K_z^u$ therefore requires the vertical gradient of the horizontal velocity (the shear), the gradient Richardson's number ($Ri$), the critical Richardson's number ($Ri_c = 0.4$) and the mixing length $l$. To ensure stability, the model sets a lower bound on the vertical shear $\lVert \partial_z {\textbf u} \rVert^2=10^{-3}$\,\unit{s^{-1}}.  Note that all these quantities can be diagnosed from the prognostic variables provided by the dynamics. 

The mixing length is modeled by 
%
\begin{equation} \label{mixing_length}
    l(z) = \left(\lambda^{-1} + \frac{1}{\kappa(z+z_0)}\right)^{-1},
\end{equation}
%
where $\kappa$ is the Von~Karman constant and $\lambda$ is the minimum mixing length. We set $\kappa=0.4$ and $\lambda = 100$\,\unit{m}.

The gradient Richardson's number~\citep{Frederic_these} is approximated as,
%
\begin{equation} \label{Richardson}
    Ri = \frac{\frac{g}{\theta}\frac{\partial \theta}{\partial z}}{\left|\left|\frac{\partial \textbf{u}}{\partial z}\right|\right|^2},
\end{equation}
%
the ratio between the buoyant production and consumption of turbulence (represented by the numerator) and the wind shear production of turbulence. $Ri_c$ is a dynamic stability criteria: when $Ri < Ri_c$, there is a net production of turbulence. The critical Richardson number is set to $Ri_c = 0.4$.  When $Ri > Ri_c$ a minimum eddy viscosity is imposed~\citep{Frederic_these},
%
\begin{equation} \label{emin_coef}
    K_z^u = l(z)\sqrt{e_{\mathrm min}}.
\end{equation}
%
where $e_{\mathrm min}$ is interpreted as a lower bound on the TKE.  The minimum TKE is set to the small value $e_{\mathrm min}=10^{-16}$\,\unit{m^2\,s^{-2}} to avoid spurious vertical diffusion in inactive regions. The maximum diffusivity $K^u_{\mathrm max}$ is set to 1~\,\unit{m^2\,s^{-1}}.

The physics model for vertical diffusion is therefore simply given by Fick's second law,
%
\begin{equation}
     \frac{\partial \textbf{u}}{\partial t} = 
     \frac{\partial}{\partial\zeta} \left( K_z^u \frac{\partial \textbf{u}}{\partial\zeta}\right), \label{turb_momentum}
\end{equation}
%
where the eddy viscosity $K_z^u$ is computed using the mixing length model (\ref{eq:eddy_viscosity}) and $\zeta$ is a floating vertical coordinate.  As mentioned above, the diffusion equation dissipates resolved kinetic energy $K$ provided by the dynamics step at the rate $\epsilon$ and diffuses this kinetic energy vertically. 

The SGS physics model for heat includes a radiative flux term and a source term $\epsilon$ due to heating (due to the dissipation of TKE and, indirectly, from the dissipation of $K$) in addition to the turbulent diffusion term,
%
\begin{equation}
     h_\theta\frac{\partial \theta}{\partial t} = 
     \frac{\partial}{\partial\zeta} \left( h_\theta K_z^\theta \frac{\partial \theta}{\partial\zeta}\right) 
     - \frac{\partial F_{rad}}{\partial \zeta} 
     + \epsilon, \label{turb_heat}
\end{equation}
%
where $\theta$ is the potential temperature, $h$ is the enthalpy, $h_\theta = \partial h/\partial \theta$.  Recall that we assume that $K_z^\theta = K_z^u$. 

The equation for resolved kinetic energy $K = 1/2 \|{\textbf u}\|^2$ is derived from the prognostic equation for velocity~(\ref{turb_momentum}) by taking its inner product with $\textbf{u}$ and rearranging,
%
\begin{equation} \label{TKE}
    \frac{\partial K}{\partial t} + \frac{1}{\rho} \frac{\partial \textbf{F}_K}{\partial \zeta} = - \epsilon
\end{equation}
%
where $\textbf{F}_K = \textbf{u} \cdot (\frac{1}{\rho} F_u \hat{{\textbf z}})$ is the vertical flux of $K$. This provides an expression for the TKE dissipated to heat, $\epsilon$, which is included in the turbulent diffusion equation for potential temperature (\ref{TKE}),
%
\begin{equation}
    \epsilon =  - (\frac{1}{\rho} F_u \hat{{\textbf z}}) \cdot \frac{\partial \textbf{u}}{\partial z} = K_z^u \left|\left|\frac{\partial \textbf{u}}{\partial z}\right|\right|^2.
\end{equation}
%
Since our turbulence model is only a first order closure approximation, there is no need for a model equation for $\epsilon$: it is diagnosed directly from the prognostic variables. It is assumed that the part of the resolved kinetic energy that is dissipated at rate $\epsilon$ is instantaneously converted to heat (rather than first being converted to TKE). The heat equation (\ref{turb_heat}) must therefore include a corresponding source term $\epsilon$. 

In contrast, in the second-order $k$-$\epsilon$ model, the eddy viscosity $K_z^u \propto k^2/\epsilon$, where $k$ is the (unresolved) small scale TKE. $K_z^u$ therefore depends on two unresolved quantities that must be modeled. They are computed by solving two coupled time-dependent model equations for their tendencies.

With a block cell vertical discretization, as explained in section Section~\ref{discrete}, the tendency equation for momentum, discretized at cell layer center $k$ is
 %
\begin{align}
     D_t\textbf{u}_k &= -\frac{1}{\rho_k} \frac{\textbf{F}_{k+1/2} - \textbf{F}_{k-1/2}}{\Delta z_k}, \\
    \textbf{F}_{k+1/2} &= \left(K^u \frac{\partial u}{\partial z}\right)_{k+1/2} = K_{k+1/2}^u \frac{u_{k+1}-u_k}{\frac{1}{2}(\Delta z_{k+1} + \Delta z_k)}\\
    \textbf{F}_{k-1/2} &= \left(K^u \frac{\partial u}{\partial z}\right)_{k-1/2} = K_{k-1/2}^u \frac{u_k-u_{k-1}}{\frac{1}{2}(\Delta z_k + \Delta z_{k-1})},
\end{align}
%
where $D_t\textbf{u}_k$ is the material derivative of momentum and $\textbf{F}_{k\pm 1/2}$ are the vertical turbulent momentum fluxes at the upper and lower interfaces $k\pm1/2$. The tendency equation for heat follows the same scheme as the momentum equation, 
 %
\begin{align}
     h_{\theta} D_t\theta_k &= -\frac{1}{\rho_k} \frac{\textbf{F}_{k+1/2} - \textbf{F}_{k-1/2}}{\Delta z_k} \\
    F_{\theta ,k+1/2} &= \left(K^u \frac{\partial \theta}{\partial z}\right)_{k+1/2} = K^u_{k+1/2} \frac{\theta_{k+1}-\theta_k}{\frac{1}{2}(\Delta z_{k+1} + \Delta z_k)}\\
    F_{\theta ,k-1/2} &= \left(K^u \frac{\partial \theta}{\partial z}\right)_{k-1/2} = K^u_{k-1/2} \frac{\theta_k-\theta_{k-1}}{\frac{1}{2}(\Delta z_k + \Delta z_{k-1})},
\end{align}
%
where the turbulent dissipation rate $\epsilon$ is neglected and $\textbf{F}$ is the total heat flux, incorporating both the turbulent and radiative heat flux. These equations are completed by Neumann (flux) boundary conditions at the bottom and top of each vertical column. The flux at the top of the atmosphere is assumed to be zero for both velocity and potential temperature, $F_{u,N} = F_{\theta,N} = 0$, while the fluxes at the ground, $F_{u,0}, F_{\theta,0}$ are provided by the surface flux parameterization scheme discussed in the following section.

The radiative heat flux is computed by the radiative transfer scheme and is an input to the turbulence model. During a time step the inputs to the model are the surface temperature and surface heat flux at the surface boundary, while the outputs of the model are the tendencies for velocity and potential temperature due to the parameterized SGS processes.

\subsubsection{Surface flux scheme}
%
The surface heat and momentum drag flux are computed using the inputs of $z, \theta$ and $u$, at the lowest (surface) layer center of the column ($k=1$). The parameterization for the fluxes are the bulk formulae of~\citet{Louis_1979}, which are based on the Monin--Obukhov similarity theory for boundary layer turbulence that describes vertical profiles of velocity and potential temperature as a function of height $z$ and the Obukhov length $L$, $z/L$. The flux model uses the following parameters, derived by \citet{Louis_1979}:

%
\begin{minipage}{0.35\linewidth}
\begin{align*}
    a &= \frac{\kappa}{\log{(z/z_0)}}\\
    b &= 4.7\\
    C^u &= 7.4\\
    C^\theta &= 5.3\\
    Ri_B &= \frac{gz(\theta - \theta_0)}{\theta (u^2+e_{min})}\\
    R &= 0.74,
\end{align*}
\end{minipage}
\begin{minipage}{0.5\linewidth}
    \begin{align*}
        c^* &= 2bC^*a^2\left(\frac{z}{z_0}\right)^{1/2}\\
        f^* &= 1 - \frac{2b Ri_B}{1 + c^*\sqrt{-Ri_B}} \tag{\textrm{if $Ri_B < 0$}}\\
        f^* &= \frac{1}{(1+bRi_B)^2} \tag{\textrm{if $Ri_B > 0$}}\\
        C_D^u &= a^2uf^u\\
    C_D^\theta &= \frac{a^2u}{R}f^\theta
    \end{align*}
\end{minipage}\\
%
\noindent where $z_0=0.1$~m is roughness length, $\theta_0$ is the potential temperature at $z_0$, $Ri_B$ is the bulk Richardson's number, $\kappa=0.4$ is the Von~Karman constant. The exponent * replaces $u$ or $\theta$ and $C_D^\theta$ are drag coefficients for velocity and potential temperature. $a$ represents the drag coefficient in the neutral surface condition case, while the coefficient calculated for $Ri_B<0$ and $Ri_B>0$ models the drag coefficients in the unstable and stable cases respectively. $R$ is an empirical parameter that relates $C_D^\theta$ to $C_D^\theta$. 

The surface fluxes are therefore computed as surface drags using the drag coefficients computed above,
\begin{align}
    \textbf{F}_u &= -\rho C_D^u \textbf{u},\\
    F_\theta &= -\rho C_D^\theta (\theta - \theta_0).
\end{align}
%

\subsubsection{Soil}\label{soilexplain}
The inclusion of a multi-layer soil model is rare in simple physics models. However, our model includes a soil column which stores and releases heat over a wide range of time scales in order to provide more accurate and realistic conditions at the surface boundary. In this model the only prognostic variable of the soil is its depth-dependent temperature $T_{s}(z,t)$, with $z<0$. The heat flux in the soil is
%
\begin{equation} \label{heatflux}
    F_s = - k \frac{\partial T_s}{\partial t}
\end{equation}
%
where $k$ is the thermal conductivity in W [K m]$^{-1}$.

From Fick's second law the temperature tendency in the soil is given by the diffusion equation~\citep{Frederic_these},
%
\begin{equation}
    \frac{\partial T_s}{\partial t} =  -\frac{1}{I} \frac{\partial F_s}{\partial z} = \frac{k}{I} \frac{\partial^2 T_s}{\partial z^2}, 
\end{equation}
%
where $k$ is the thermal conductivity and the thermal inertia $I = 3000$\,\unit{J\,m^2\,K^{-1}\,s^{-1/2}}.

The soil is discretized by $N_S=10$ layers, where the heat fluxes are computed at layer interfaces $l=-N_S,\ldots,0$, where $l=0$ is the surface. The interfaces are located at generalized depths 
%
\begin{align}
    \zeta'_l = \frac{\alpha^{l} - 1}{\alpha - 1}\sqrt{T},
\end{align}
%
where the time scale $T = 2000$\,\unit{s}$/\pi$ and exponent $\alpha=2$ are model parameters. Depth is therefore measured as in terms of time, where deep layers correspond to increasingly slow heat diffusion time scales. Including more vertical layers add more slow time scales to heat diffusion in the soil column.  Table~\ref{tab:layer_depth} shows the time scales corresponding to each soil layer for the $N_S=10$ soil layer case considered in this paper.
%
\begin{table}[h!]
\centering
\begin{tabular}{|c|r|r|}
\hline
\textbf{Layer} & \textbf{Period}     & \textbf{Depth [m]} \\
\hline
1  & 5.6 hours    & 0.11  \\
2  & 2.1 days    & 0.33  \\
3  & 11 days    & 0.76  \\
4  & 52 days    & 1.6   \\
5  & 220 days    & 3.4   \\
6  & 2.5 years     & 6.9   \\
7  & 10 years    & 14  \\
8  & 41 years    & 28  \\
9  & 170 years & 56  \\
10 & 660 years & 110 \\
\hline
\end{tabular}
\medskip
\caption{The time periods of each soil layer and their equivalent depth.  10 soil layers is sufficient for climate simulations of up to about 1000 years. The calculated depth in the table uses the heat capacity $C_p=2.2 \times 10^6$\,\unit{J\,m^{-2}\,K^{-1}} and thermal conductivity $\lambda=4.09$\,\unit{J\,s^{-1}\,m^{-2}\,K^{-1}}.}
\label{tab:layer_depth}
\end{table}

The first interface below the surface is ($\zeta'_1 = \sqrt{T}$) and it includes the damping depth of the diurnal wave. The damping depth is the depth of a wave of period $P$ whose amplitude is reduced by $e$~\citep{IntroMicrometeorology}. More generally, the depth of penetration of a wave of period $P$, in days, is $\sqrt{P}$~\citep{Frederic_these}. More layers included in the soil model allows for the release of energy at more (longer) time scales, for example a diurnal, monthly and annual time scales. 

Fluxes at interface $l$ are computed by:
%
\begin{align}
    F_{s,l} = I \frac{T_{s,l+1/2} - T_{s,l-1/2}}{\frac{1}{2}(\Delta \zeta'_{l+1/2} + \Delta \zeta'_{l-1/2})},
\end{align}
%
where $\Delta \zeta'_{l+1/2} = \zeta'_{l+1} - \zeta'_{l}$. Therefore, the soil temperature tendency of soil layer $l+1/2$ is defined by
%
\begin{align}
    \frac{\partial T_s}{\partial t} = -\frac{1}{I} \frac{F_{s,l} - F_{s,l-1}}{\zeta'_l - \zeta'_{l-1}}.
\end{align}
%
If there are $L_s$ layers, the Neumann boundary conditions for the soil column are
%
\begin{align}
    F_{s,0} & = F_{r} + F_H, \tag{\textrm{Flux at ground boundary, $l=0$}}\\
    F_{s,L_s} &= 0, \tag{\textrm{Flux of bottom soil interface} $l=L_s$}
\end{align}
%
where $F_{r}$ is the surface radiative flux and $F_H$ is the surface heat flux.

\subsubsection{Dry Convective Adjustment}
%
Dry convection must be parameterized because the dynamics makes the hydrostatic approximation. The vertical turbulence model only simulates turbulent diffusion, not convection. It could therefore produce an unstable vertical temperature stratification, characterized by $\partial \theta/\partial z < 0$~\citep{Frederic_these}. The dry convective adjustment scheme of the physics model mitigates any instability by instantaneously relaxing the temperature to an adiabatic profile. 

A layer is unstable if its potential temperature is greater than that of the layers above.  The convection parameterization scheme checks each column to see if a layer is unstable. If a particular layer of the one-dimensional vertical column is unstable, local mixing is applied. This local mixing is modeled by a temperature relaxation to the mass weighted temperature average of the unstable layers. This scheme also takes into account the transport of momentum that occurs with convection by mixing the momentum in the unstable layers. For each unstable layer in a column only a proportion $\alpha$ of the cell is mixed, where
%
\begin{equation}
    \alpha = \frac{\int |\overline{\theta} - \theta| \,dp}{\int \theta \,dp},
\end{equation}
%
and $\alpha < 0$, which is always verified \citep{Frederic_these}.


\subsection{Physics package usage, input and output}\label{Package_Inputs}
The vertical physics model is intended to be coupled with a dynamical core. This coupling is usually done via a split step method, where the ``dynamics'' and ``physics'' components take alternate times steps. The prognostic variables of the physics package are temperature, zonal velocity and meridional velocity. When a physics step is taken, the package outputs the tendencies of the prognostic variables.  Note that each of the the physics sub-models (radiation, turbulent diffusion, dry convective adjustment, diurnal cycle) can be turned on or off separately.

One of the major objectives of the package is its versatility in modeling different planets. In order to accommodate this goal, the package initialization allows users to read in input parameters to tailor the simulation to the desired planet. Beyond these parameters, calls to initialization routines set up the logistics of the grid. 
%
\subsection {Numerical integration and boundary coupling}
%
\subsubsection{Time integration\label{sec:physics_time}}
%
To solve the system of differential equations at the next time step a numerical integration method is employed. The explicit Euler method is a first-order method commonly used to numerically integrate a system of differential equations. Given an autonomous system of ODEs,
%
\begin{equation}
    \frac{d\mathbf{S}}{dt}=\mathbf{F}(\mathbf{S}(t)),
\end{equation}
%
where $\mathbf{F}(\mathbf{S}(t))$ is called the ``tendency''. The stencil of the Explicit Euler method is given by the linear Taylor series approximation of $\mathbf{S}(t*) = \mathbf{S}(t+\tau) + O(\tau^2)$
%
\begin{equation}
    \mathbf{S}(t*) = \mathbf{S}_t + \tau \mathbf{F}(\mathbf{S}_t)
\end{equation}
%
where:
%
\begin{itemize}
    \item $t$ is the current time
    \item $\tau$ is the time step (assumed fixed)
    \item $t*$ = t + $\tau$ is the time at the next discrete time
    \item $\mathbf{S}_t$ is the numerical approximation of $\mathbf{S}(t)$, $\mathbf{S}_t\approx\mathbf{S}(t)$.
\end{itemize}

A variety of time integration methods have been developed from the Explicit Euler method by making different approximations for the tendency of $\mathbf{S}$, $\mathbf{F}(\mathbf{S}$),  (e.g.\ Runge-Kutta schemes). These methods differ in their order of accuracy and numerical stability properties. In the Explicit Euler method the tendency is simply approximated as $\mathbf{F}(\mathbf{S}_t)$, where the values of $\mathbf{S}$ are at time t. The implicit Euler method improves numerical stability by using the tendency for \textbf{S} at the the next time step ($\mathbf{F}(\mathbf{S}_{t*})$). This requires more computation since a linear system must be solved at each time step, in general. However, Implicit Euler has a much wider stable domain compared to the Explicit Euler method, allowing for a larger time step and increasing its desirability. The radiation model uses Explicit Euler, while the soil and turbulence models use Implicit Euler. 

%The variation of Euler's method that is used in this physics package is semi-implicit Euler's method where the tendency is $\mathbf{F}(\mathbf{S}_{t}, \mathbf{S}_{t*})$. This is more complex as the tendency is a function of both $\mathbf{S}_t$ and $\mathbf{S}_{t*}$.

\subsubsection{Boundary coupling}
%
Since the sub-components of the physics are air and soil, it is essential to specify how they are coupled in the time integration scheme. The air and soil columns are coupled via their boundaries, which is the horizontal interface between the soil and air (i.e., the ground), and the boundary conditions need to be taken into consideration and coupled when evaluating the tendencies of the two components. The simplest boundary coupling is explicit, and the boundary conditions  are evaluated at time $t$ from the previous state. Therefore, the time integration equation includes the boundary conditions $\mathbf{q_t}$,
%
\begin{equation}
    \mathbf{S}_{t*} = \mathbf{S}_t + \tau \mathbf{F}(\mathbf{S}_t,\mathbf{q}_t).
\end{equation}
%
Since the air and soil are strongly coupled, coupling the boundary conditions implicitly is desirable for improved stability. (Note that both the soil and air physics strongly affect temperature.)  Therefore, the boundary conditions are evaluated at time $t*$. The tendency is approximated as 
%
\begin{equation}
    \mathbf{F}(\mathbf{S}(t),\mathbf{q}(t)) \approx 
    \mathbf{F}(\mathbf{S}_t,\mathbf{S}_{t*},\mathbf{q}_{t},\mathbf{q}_{t*})
\end{equation}
%
\begin{equation}
    \mathbf{S}_{t*} = \mathbf{S}_t + \tau \mathbf{F}(\mathbf{S}_t,\mathbf{q}_t) + \frac{\partial \mathbf{S}_{t*}}{\partial \mathbf{q}_{t*}} (\mathbf{q}_{t*} - \mathbf{q}_t),
\end{equation}
%
where $\partial \mathbf{S}_{t*}/\partial \mathbf{q}_{t*}$ represents the dependence of $\mathbf{S}_{t*}$ on the boundary conditions at the new time $t^*$. 

\subsection{Spatial discretization}\label{discrete}
%
The physics package models the effects of the small scale unresolved physics on each 1D vertical column on the sphere, taking the current state of the prognostic dynamical variables as inputs. The vertical discretization is a block cell discretization, with all prognostic variables located at the center of each layer cell, while the fluxes are located at the interfaces between each layer. Each column has $N$ layers, which is an input parameter to the package from the dynamics. There are a total of $N + 1$ interfaces, where the surface (ground) interface has index 0 and top interface (top of atmosphere) has index $N$. The soil column is considered separately, with negative indices. With the chosen floating vertical coordinates, the height of each cell layer ($\Delta z_k$) changes over time. Technical details are provided in~\ref{sec:software}.

\subsection{Limitations, Decisions and Possible Extensions} \label{limits_decisions}
%
The simple dry physics package makes approximations to ensure it is relatively simple, stable and easy to use and understand. However, these choices necessarily mean the physics is idealized and not sufficiently realistic for many applications.

First, a dry climate is assumed by the package which neglects condensible substances (e.g., water or methane) on the atmosphere. The lack of moisture is a major limitation, especially for planets with condensible substances that have a major effect on the atmosphere, like Earth. However, assuming a dry climate removes the added complexity of the effects of condensible substances and makes the physics much simpler. Due to this limitation, we cannot expect accurate results when modeling Earth. For example, the temperature of the  surface temperature should be much higher than that of Earth. This is due to the fact that the lapse rate of the simulation is constrained to be close to the dry adiabatic lapse rate, to keep a physically consistent dry model, while the lapse rate for Earth is closer to the moist lapse rate.

Secondly, the surface flux parameterization in the physics package was created empirically using Earth observation data, which includes moisture. Therefore the scheme may not be suitable for modeling dry planets and might require tuning when modeling planets whose condensible substances are not water.

Finally, the current radiation parameterization uses only two frequency bands: shortwave (visible) and longwave (infrared). In general, radiation models are  expensive, and therefore having only two frequency bands keeps the cpu time per step minimal. However, with only two radiation bands, when modeling Earth, the temperature will not be accurate in the stratosphere as UV light has a major effect on the temperature.   Nevertheless, unlike for moisture, it is relatively simple to include additional radiation bands.

Adding a third absorption band for UV light would improve the realism of the vertical temperature profile in the stratosphere, producing the expected trend of an increase in temperature with height. However, this addition would increase the computational cost of the radiation step.

\section{Configuration and Numerical Implementation\label{Configuration}}
The simple dry physics package needs to be coupled with a dynamical core to simulate the atmospheric general circulation of a planet. In this section we describe the {\sc wavetrisk} dynamical core~\cite{wavetrisk} used to evaluate the characteristics of the physics package, the configuration of the physics-dynamics coupling and its numerical implementation. We also present the base simulation of the~\cite{HeldSuarez} dry physics model, which is to be used for future comparison of the simple dry physics package.

\subsection{{\sc wavetrisk} dynamical core}
Since we want to better understand the need for scale-aware physics, and how dynamic adaptivity combines with a physics package, we use the {\sc wavetrisk} dynamical core. {\sc wavetrisk} is a global three-dimensional dynamically adaptive hydrostatic dynamical core that can be run in both adaptive and non-adaptive modes (as well as in atmosphere or ocean configurations). This allows us to directly compare equivalent adaptive and non-adaptive runs. {\sc wavetrisk} utilizes a hexagonal-triangular C-grid with the TRiSK horizontal discretization scheme, which is well known for its mimetic properties and stability.  

Dynamic horizontal grid adaptivity is implemented using a multiscale wavelet approach with user-specified minimum (coarsest) and maximum (finest) resolutions. Nodes and edges are added or removed at each scale as necessary to achieve a user-specified relative error tolerance for each prognostic variable (i.e., pseudo-density $\mu = \rho_0 \Delta z$, mass-weighted potential temperature $\Phi = \mu\theta$, and edge velocities $u,v,w$). A node/edge is retained if it is retained in any vertical layer.  The computational grid is adapted after each full time step as follows: 
\begin{enumerate}
\item 
Compute wavelet coefficients at each horizontal position/vertical layer/horizontal resolution scale for all prognostic atmosphere variables.
\item 
If a wavelet coefficient is above the dimensional tolerance threshold for any prognostic variable at any vertical layer, keep the associated grid point. If not, delete it. 
\item 
Add to the active grid the nearest neighbour grid points at the same scale and at the finer scale. This allows for a local scale-dependent CFL criterion $u\Delta t/\Delta x\leq 1$ and accounts for quadratic nonlinearities (which can generate scales twice as small in a single time step).
\item 
Interpolate all prognostic variables (including soil temperature) onto the new adapted grid. 
\end{enumerate}
This process ensures that the grid adaptation algorithm conserves mass and various mimetic properties. The resulting adaptive grid is a set of columns of varying cross-section. {\sc wavetrisk} may be run either adaptively with a range of horizontal resolutions, or non-adaptively with a fixed horizontal grid resolution. 
Full details on the {\sc wavetrisk} model are given in~\cite{Kevlahan/Dubos:2019}. 

\subsection{Model configurations}
Time stepping for the physics-dynamics coupled model uses explicit low storage third-order Runge-Kutta scheme for {\sc wavetrisk} and an implicit Euler step for the physics package. Both the dynamical core and the physics package have the same time step, with {\sc wavetrisk} determining the time step and supplying it to the physics. The adaptive simulation uses the time step determined by the finest resolution.

We use horizontal bi-Laplacian hyperdiffusion for vorticity $\nabla\times\mathbf{u}$, velocity divergence $\nabla\cdot\mathbf{u}$ and scalars to reduce grid scale noise. The hyperviscosity $\nu_4$ is defined non-dimensionally as a proportion of the maximum stable value for the given time step $\Delta t$ and minimum average grid size $\Delta x_{\mathrm{min}}$ ($\Delta t$ is based on $\Delta x_{\mathrm{min}}$).  The maximum stable viscosity is the analytical value for the hexagonal/pentagonal C-grid on the plane, modified with empirical correction factors for the irregular multiscale hexagonal C-grid on the sphere. (The empirical corrections are conservative by 15\% to 33\%, and are based on Gershgorin circle theorem estimates evaluated on the actual multiscale grids.) Note that the TRisK discretization~\citep{Ringler/etal:2010} is stable without horizontal diffusion on a fixed resolution grid, but there is some grid scale noise.

{\sc wavetrisk} uses Lagrangian (floating) vertical coordinates with periodic conservative remapping of the prognostic variables to the initial hybrid pressure grid. Note that vertical remapping introduces some vertical diffusion. The data structure combines the atmosphere layers and soil layers into a single vertical column, but only the atmosphere layers require remapping. All simulations have 30 hybrid pressure layers, with the top of the atmosphere at 225~Pa.  The $A,B$ hybrid pressure coefficients are given in~\cite{Ullrich2012_DCMIP_TestCases} Table~XVIII. Layer~6 is at 850~hPa and layer 12 is at 533~hPa.  There is no orography. There are 10 soil layers for all simulations.



\begin{table*}[t]
\caption{Average horizontal grid resolutions $\Delta~x$in degrees and km  in kilometers.}
\label{ResolutionRange}
\begin{tabular}{lcr}
\tophline
Resolution & $\Delta x$ (km)\\
\middlehline
$2^\circ$ & 240\,\unit{km}\\
$1^\circ$ & 120\,\unit{km}\\
$0.5^\circ$ & 60\,\unit{km}\\
\bottomhline
\end{tabular}
\belowtable{} % Table Footnotes
\end{table*}

The physics uses Earth astronomical values and boundary drag parameterizations, including seasons (see Table~\ref{tab:astronomical}). The initial conditions for prognostic variables are based on an isothermal temperature profile of $250$\,\unit{K} and zero velocity. While it would be more realistic to have an initial temperature that depends on latitude, a horizontally uniform isothermal profile is a simple initial condition that has been used in Dynamical Core Intercomparison Projects (DCMIP) (e.g. \citet{DCMIP_ic}) and the climate dynamics develops within a few days.

The reference resolution used to characterize the climatology of the {\sc wavetrisk}--simple physics coupling is a relatively coarse non-adaptive $2^\circ$ ($\sim$~240\,\unit{km}) simulation. A resolution of $2^\circ$ is considered sufficient for coarse climate simulations.

The adaptive simulation uses a minimum resolution of $2^\circ$  ($\sim$~240\,\unit{km}) and a maximum resolution of $ 0.5^\circ$ (60\,\unit{km}). This means that the maximum possible refinement is 4 times and the maximum grid compression is $4^2 = 16$ (i.e., about 16 times faster than a non-adaptive simulation with the maximum resolution).  The simulations are initialized at the coarsest $2^\circ$ resolution and refined and coarsened after each time step with a relative error tolerance of 5\% for based on the $\ell_2$ norms for each vertical layers.

\subsection{Simulation and diagnostic configuration}
We first compare a non-adaptive $2^\circ$ {\sc wavetrisk}--simple physics simulation without seasons to an equivalent Held \& Suarez~\citep{HeldSuarez} simulation and then include seasonal radiation forcing (all simulations include diurnal radiation forcing).  The non-adaptive $2^\circ$ climatology is the reference for the higher resolution $0.5^\circ$ non-adaptive and adaptive simulations. 

All simulations are run for five (Earth) years. This period includes two years for the means of the prognostic variables in all layers to converge to within about 2\% for runs with seasons.  The value of 2\% was chosen (rather than a smaller value) since the climate variables have significant intrinsic variance due to turbulence. %After three years all layers are within 0.2\% of their means.

The climatology of the physics-dynamics coupled simulations is characterized by 5-year zonal first- and second-order statistics: mean temperature, mean zonal velocity, eddy kinetic energy, eddy heat flux temperature variance and eddy momentum flux. In addition, we include the vertical profile of temperature, averaged over the sphere and instantaneous horizontal slices of temperature and vorticity. The results for the $0.5^\circ$ adaptive simulation include a visualization of instantaneous three-dimensional vorticity isosurfaces, zonally averaged vertical velocity and energy spectrum results.  

\subsubsection{Held-Suarez base simulation} \label{HS-base}
The \ac{HS} model \citep{HeldSuarez} is the simplest dry physics model that can be coupled to dynamical cores and a standard benchmark used for comparison of other simplified physics models. Therefore our initial comparison for the simple dry physics package will be to this model which includes only a relaxation of temperature to a radiative equilibrium value and a layer dependent linear damping of velocities in the surface boundary layer.

The \ac{HS} simulation used for future comparison is from a non-adaptive {\sc wavetrisk} simulation with maximum resolution $1^\circ$, using the initial conditions of \citealt{Jablonowski2006}. Figure \ref{HS_climatology} presents the zonal statistics of the 5-year simulation. (Note that~\cite{wavetrisk} did not find significant qualitative or quantitative differences in the climatology when maximum horizontal resolution is increased to $1/4^\circ$.)

\begin{figure*}[t]
\includegraphics[width=\linewidth]{Figures/HSJ6Z30_climatology.png}
%
\caption{Climatology from a {\sc wavetrisk}--Held \& Suarez physics simulation with uniform $1^\circ$ resolution and 30 vertical layers. Statistics are computed from days 200 to 1200. The first order zonal statistics presented are the (a) temperature and (b) zonal velocity. The second order zonal statistics are (c) eddy kinetic energy, (d) eddy heat flux, (e) temperature variance and (f) eddy momentum flux.}
\label{HS_climatology}
\end{figure*}

It is important to note that, although it is a standard benchmark, and therefore an interesting comparison, the \ac{HS} model is fundamentally different from the simple physics model, and so we should not expect qualitatively or quantitatively similar statistics.  In the \ac{HS} model temperature is relaxed exponentially to a specified latitude and height-dependent profile, based on Earth climatology, and Rayleigh friction damps velocity near the surface (i.e., for pressures greater than 700\,\unit{hPa}) with a 1~day relaxation time.  There are no other physical processes represented and the physics does not directly modify the dynamics, except by through the temperature relaxation and boundary friction.  In addition, the simple physics model includes diurnal radiation forcing and there is no diurnal scale in the \ac{HS} temperature relaxation (the relaxation time scale is about 40~days in the stratosphere). 

\section{Results\label{Results}}
\subsection{Non-adaptive $2^\circ$ climatology without seasons}\label{NA-climatol-no-seasons}
%
Figure~\ref{TempProfileClimatologyNonAdapt} presents the temperature profile for the non-adaptive simple physics simulation without seasons compared to the \ac{ISA} model profile and the equilibrium state of the \ac{HS} model~\citep{HeldSuarez}, explained in the previous section. The \ac{ISA} is a standard temperature profile for Earth. 

The simple physics profile, plotted in red, is spatially averaged over the sphere at each layer and temporally averaged over five years. The temperature profile decreases monotonically with height, without a well-defined tropopause. In comparison to the standard atmosphere and \ac{HS}, the simple physics temperature profile has a significantly warmer surface and cooler upper stratosphere. A hotter surface is expected to the lack of moisture in the atmosphere, while the continued decrease in upper atmosphere temperature is due to the lack of a UV radiation absorption band in the stratosphere, as mentioned in section \ref{limits_decisions}. 
%
\begin{figure*}[t]
\centering
\includegraphics[width=0.6\textwidth]{Figures/TempProfileClimatology_J5Z30.pdf}
\caption[Climatology Vertical Profile of Temperature Averaged over five years]{Temperature profile averaged over five years from a coupled {\sc wavetrisk} and simple physics simulation with 10 soil layers. The observed trend of the profile is a decrease in temperature with height. The figure also displays the \ac{ISA}  and a $1^\circ$ 30 layer Held--Suarez simulation profile for comparison.}
\label{TempProfileClimatologyNonAdapt}
\end{figure*}

It is helpful to begin by comparing the simple physics climatology without seasons to the most basic physics model of ~\cite{HeldSuarez}. To do this, we ran a five year simple physics coupled to {\sc wavetrisk} non-adaptive simulation without seasons (i.e., obliquity of $0^\circ$). All other parameters were as described above and standard zonal statistics are presented. The {\sc wavetrisk}--\ac{HS} simulation is as presented in section \ref{HS-base}. When compared to the \ac{HS} results shown in Figure~\ref{HS_climatology}, the simple physics climatology in Figure~\ref{ZonalStatsNoSeasonNonAdapt} has significant differences, although they are broadly similar. The biggest differences are in the eddy \ac{KE}, which have a different zonal jet shape, and in the zonal velocity which has its highest velocity near the top of the atmosphere, rather than at a similar pressure as the highest eddy \ac{KE}. Furthermore, while the eddy heat flux and temperature variance differs between the simple physics and the \ac{HS} simulations, the overall structure of the two statistics are similar, with high values near the mid-latitudes and at a normalized pressure of $\sim0.85$. 

\begin{figure*}[t]
\includegraphics[width=\linewidth]{Figures/SimpleJ5Z30_climatology_noseasons.png}
\caption{Climatology of the zonal statistics from a {\sc wavetrisk}--simple physics (non-adaptive) simulation without seasons. Horizontal resolution is $2^\circ$, there are 10 soil layers and 30 vertical atmosphere layers. The first order zonal statistics presented are the (a) temperature and (b) zonal velocity. The second order zonal statistics include (c) eddy kinetic energy, (d) eddy heat flux, (e) temperature variance and (f) eddy momentum flux. All statistics were time averaged over 5 years.}
\label{ZonalStatsNoSeasonNonAdapt}
\end{figure*}

\subsection{Non-adaptive $2^\circ$ simulation climatology including seasons}\label{NA-climatol-with-seasons}
%
We now re-run the previous non-adaptive simulation, but this time including seasons (i.e., obliquity $23.5^\circ$). Figures~\ref{ZonalStatsNonAdapt_SpringFall} and \ref{ZonalStatsNonAdapt_SummerWinter} presents the vertical profile of the five year seasonal climatology for the first- and second-order zonal statistics. Unsurprisingly, the statistics are markedly different with seasons compared to the equivalent results without seasons shown in~\ref{ZonalStatsNoSeasonNonAdapt}. As pressure decreases with height, the temperature decreases, but there is also minimal zonal temperature variance near the top of the atmosphere. The large temperature variance at the poles is simply due to the effect of seasons. Furthermore, the eddy \ac{KE} presents a profile that is indicative of a zonal jet at the mid latitudes with a peak variance at a normalized pressure of $\sim~0.25$.  Overall, simple physics with seasons has qualitatively and quantitatively distinct climatology from Held--Suarez, with only the mean zonal temperature and high altitude peaks of eddy kinetic energy at $\pm 40$~N/S and $P/P_S = 0.25$ qualitatively similar.

An interesting feature of the simple physics with seasons climatology is the peaks in eddy \ac{KE}, eddy momentum and heat flux at $15^\circ$~N/S near the surface. This feature does not appear in the simple physics simulation without seasons shown in Figure~\ref{ZonalStatsNoSeasonNonAdapt}. This suggests that seasonal radiation forcing greatly changes the climatology, both (unsurprisingly) at the poles, but also in the surface boundary layer near the equator. The Held--Suarez model cannot capture this effect.
%
\begin{figure*}[t]
\centering
{\bf Summer}\\
\includegraphics[width=\textwidth]{Figures/SimpleJ5Z30_climatology_summer.png}
{\bf Fall}\\
 \includegraphics[width=\textwidth]{Figures/SimpleJ5Z30_climatology_fall.png}
\caption{Summer and and fall climatology of zonal statistics from a coupled $2^\circ$ resolution {\sc wavetrisk}--simple physics (non-adaptive) simulation with seasons ($23^{\circ}$ obliquity). Note much smaller temperature variance than in spring/fall (different axis limits).}
\label{ZonalStatsNonAdapt_SummerFall}
\end{figure*}

The longitude-latitude horizontal projection of five year average surface temperature climatology is shown in Figure~\ref{ClimatologyTempNonAdapt}~(top). The expected structure of a warmer equator and colder poles is confirmed. However, there is still some evidence of the impact of turbulence and seasonality in the lower temperature near $\sim 100^{\circ}$ longitude.  A typical instantaneous temperature is shown in  Figure~\ref{ClimatologyTempNonAdapt}~(middle) on March 22 of year 5, which reveals the the influence of turbulence and  east-west and north-south asymmetries in the flow (i.e., prevailing mid-latitude westerly winds advecting temperature filaments).   The turbulence structure underlying the temperature advection is revealed by the vorticity field in Figure~\ref{ClimatologyTempNonAdapt}~(bottom).  The turbulence  is characterized by intense small scale vorticity filaments in the northern hemisphere and weaker larger scale vorticity filaments in the southern hemisphere.   Comparing instantaneous temperature and vorticity projections, it is clear that the temperature fluctuations follow the vorticity, confirming that  temperature is behaving approximately like a passive scalar and the turbulence structure determines the temperature structure. Overall the instantaneous projections reveal the season-dependent ``weather" of the coupled climate model, which is not observed in the uncoupled runs of the dynamical core. 
%
\begin{figure*}[t]
\centering
\includegraphics[width=0.7\textwidth]{Figures/J5Z30_Temp_Vort_LonLat.png}
\caption{Longitude-latitude projections of surface layer temperature and vorticity from the non-adaptive {\sc wavetrisk}--simple physics simulation with seasons. Top: five year mean temperature. Middle: temperature on March~22 of year~5. Bottom: vorticity on March~22 of year~5, showing turbulent structure at the surface.}
\label{ClimatologyTempNonAdapt}
\end{figure*}

\subsection{Adaptive high resolution $0.5^\circ$ climatology including seasons}
%
Dynamically adaptive GCMs have the potential to significantly improve accuracy and computational efficiency by adjusting the local horizontal grid resolution at each time step to ensure a given tolerance.  In this section we investigate our key questions, presented in section \ref{intro} regarding the dynamically adaptive \ac{GCM}, focusing on the adaptive physics-dynamics coupling's stability, efficiency to captures small-scale structure and power-law scaling, the need for adapting on the physics and restarting from a higher-resolution checkpoint. The section presents the results of the simple physics package coupled to {\sc wavetrisk-atmosphere} run adaptively with a coarsest average grid resolution of 2$^\circ$ and a finest average grid resolution of 0.5$^\circ$ (i.e., local grid refinement of up to 16 times). As for the non-adaptive case, there are 30 vertical atmosphere layers and 10 soil layers. The relative error tolerance is $5\%$ for all prognostic dynamical variables with the $\ell_2$ tolerance normalization determined separately for each atmosphere layer. The grid is adapted explicitly only on the prognostic dynamical variables, not on the physics (e.g., tendencies of the vertical diffusion or radiation sub-models). The simulation is run for five years and statistics are computed over the full five years.  

The grid of the adaptive simulation has an average compression ratio of $\sim 2.3$, meaning that there are 2.3 times fewer active cells than the equivalent finest non-adaptive grid ($0.5^{\circ}$). Over a simulation year the compression ratio varies depending on the season, with higher compression at the solstices and a lower compression at the equinoxes due to the changing weather of the season. The lowest and highest compression ratios of the simulation are about 1.8 and 3.0 respectively. Figure~\ref{fig:Adaptive_grid} displays the adaptive grids with the lowest and highest compression ratio during the spring equinox on March 22 of simulation year 5 and during the summer solstice on June 25  of simulation year 4 respectively.  At the summer solstice the northern (hotter) hemisphere is less active than the southern (cooler) hemisphere. Overall, the atmosphere is much less turbulent at the solstices than at the equinoxes. The annual pattern of grid compression is shown in Figure~\ref{fig:compression}. The atmosphere is less turbulent (higher compression ratio) at the solstices than at the equinoxes. The average compression rate for this adaptive simulation is about 2.3, and it varies between about 1.7 and 3.6.  The average compression ratio increases with the number of resolution levels.

\begin{figure*}[t]
\centering
\includegraphics[width=0.6\textwidth]{Figures/AdaptiveGrid_J5J7_T365.png}\\
\includegraphics[width=0.6\textwidth]{Figures/AdaptiveGrid_J5J7_T311.png}
\caption{Top: adaptive horizontal dual triangular grid with relative tolerance $\varepsilon = 2.5\%$ on 5 March 22 of year 5 during an equinox. The adaptive grid has a low compression ratio:  1.8~times fewer cells compared to the non-adaptive $0.5^\circ$ degree grid. Bottom: adaptive horizontal grid on June 25 of year 4 during a solstice when the adaptive grid has a high compression ratio of 3.0. At the summer solstice that the northern (hotter) hemisphere is less active than the southern (cooler) hemisphere.}
\label{fig:Adaptive_grid}
\end{figure*}

\begin{figure*}[t]
\centering
\includegraphics[width=0.8\textwidth]{Figures/compression.pdf}
\caption{Compression ratio for the adaptive simulation showing seasonal variation over 5 years. The solstices and equinoxes are labelled. The highest compression ratios are at the solstices and the lowest compression ratios are at the equinoxes.}
\label{fig:compression}
\end{figure*}
%
\begin{figure*}[t]
\centering
{\bf Summer: adaptive resolutions $2^\circ, 1^\circ, 0.5^\circ$}\\
\includegraphics[width=\textwidth]{Figures/SimpleJ5J7Z30_climatology_summer.png}
{\bf Summer: fixed resolution $0.5^\circ$}\\
\includegraphics[width=\textwidth]{Figures/SimpleJ7Z30_climatology_summer.png}
\caption{Comparison of summer climatology from the adaptive {\sc wavetrisk} simulation with relative tolerance $\varepsilon = 2.5\%$ and a non-adaptive simulation coupled to simple physics with seasons with resolutions of $2^\circ, 1^\circ, 0.5^\circ$. Statistics were averaged data saved every 5 days.}
\label{fig:ZonalStats_Summer}
\end{figure*}
%
\begin{figure*}[t]
\centering
{\bf Fall: adaptive resolutions $2^\circ, 1^\circ, 0.5^\circ$}\\
\includegraphics[width=\textwidth]{Figures/SimpleJ5J7Z30_climatology_fall.png}
{\bf Fall: fixed resolution $0.5^\circ$ }\\
 \includegraphics[width=\textwidth]{Figures/SimpleJ7Z30_climatology_fall.png}
\caption{Comparison of fall climatology from an adaptive {\sc wavetrisk} simulation with relative tolerance $\varepsilon = 2.5\%$ and a non-adaptive simulation coupled to simple physics. Statistics were averaged data saved every 5 days.}
\label{fig:ZonalStats_Fall}
\end{figure*}
%
Figures~\ref{fig:ZonalStats_Summer} and \ref{fig:ZonalStats_Fall} present the first- and second-order zonal statistics of the adaptive simulation. The statistics of the adaptive and non-adaptive simulations are qualitatively similar. The main quantitative difference is that the adaptive high resolution simulation has a more energetic eddy kinetic energy in the zonal jet structure, especially at the centre of the jet (at height $P/P_s \approx 0.25$).  The relative $l_2$ difference of the adaptive simulation compared to the non-adaptive simulation are summarized in Table~\ref{tab:adapt_errors}.  The differences are 10--18\% of the maximum values, except for zonally averaged temperature which has a much smaller error of about 0.5\%.  Figure~\ref{fig:ZonalStats_diff_SummerFall} shows that the differences are greatest at high latitudes, where the adaptive grid is coarsest around the solstices, and at high altitudes where pressure is very low (see Figure~\ref{fig:Adaptive_grid}).
%
\begin{table}[h]
\begin{tabular}{llll}
\hline
Field                & Relative $l_2$ difference (summer) & Relative $l_2$ difference (fall) & Relative $l_2$ difference (full year) \\ \hline
Temperature          & 0.50\%                        & 0.36\%                       & 0.57\% \\ 
Zonal velocity       & 15\%                          & 15\%                        & 12.8\% \\ 
Eddy heat flux       & 9.8\%                         & 14\%                        & 18\% \\
Eddy kinetic energy  & 15\%                          & 18\%                        & 11\% \\
Eddy momentum flux   & 18\%                          & 17\%                        & 9.8\% \\
Temperature variance & 12\%                          & 5.5\%                       & 0.32\% \\ \hline
\end{tabular}
\smallskip
\caption{\label{tab:adapt_errors} Relative $l_2$ zonal climatology differences between the adaptive and non-adaptive simple physics simulations. The corresponding zonal figures for summer and fall are shown in Figure~\ref{fig:ZonalStats_diff_SummerFall}. }
\end{table}

\begin{figure*}[t]
\centering
{\bf Summer}\\
\includegraphics[width=\textwidth]{Figures/Simple_climatology_diff_summer.png}
{\bf Fall}\\
\includegraphics[width=\textwidth]{Figures/Simple_climatology_diff_fall.png}
\caption{Difference between adaptive and non-adaptive summer and fall climatology of the zonal statistics. Statistics are averaged data from saved every 5 days. Note difference in scales compared to the climatology shown in Figures~\ref{fig:ZonalStats_Summer} and \ref{fig:ZonalStats_Fall}.}
\label{fig:ZonalStats_diff_SummerFall}
\end{figure*}

The adaptive simulation 5 year temperature climate normal latitude--longitude projection, seen in Figure~\ref{fig:ClimatologyTempVortAdapt} (top), displays the expected meridional temperature trends and both the climate normal and instantaneous temperature projections are qualitatively similar to the non-adaptive simulations. With the instantaneous projection taken at the end of the 5 year simulation, March 22, it displays the changing seasonal dynamics with a cooler northern hemisphere and warmer southern hemisphere. The instantaneous temperature and vorticity projections (Figure~\ref{fig:ClimatologyTempVortAdapt} middle, bottom) capture the fine scale structure, which is desired with an increase in resolution. Especially at the mid-latitudes, the eddies are more distinct in the adaptive case, thus confirming the ability of the adaptive \ac{GCM} to capture the emergence of small scale structure. 

Figure \ref{fig:TempProfileClimatologyJ5J7} compares the the vertical temperature profiles averaged over the 5 year simulation for non-adaptive simulations with resolutions of $2^\circ$ (J5) , $1^\circ$ (J6), and $0.5^\circ$ (J7) and the adaptive resolution which is a combination of the three resolutions (J5, J6, J7). All profiles present the trend of a decrease in temperature with height and among the non-adaptive (single) resolutions there is minimal difference between the profiles. However, near the top of the atmosphere, seen in Figure~\ref{fig:TempProfileClimatologyJ5J7} (b) the finest resolution ($0.5^\circ$, J7) is quantitatively closer to the the adaptive (J5-J7) profile, compared to the lower half of the atmosphere where the J7 profile is closer to the other non-adaptive resolution profiles. As mentioned earlier, it is important to note, that near the top of the atmosphere the profile does not depict Earth's trend due to the simplicity of the radiation sub-model (see section \ref{limits_decisions}). Furthermore, while the adaptive simulation profile is quantitatively near that of the three non-adaptive profiles, it is apparent that its vertical profile is cooler. 
%
\begin{figure}
    \centering
    \includegraphics[width=0.8\linewidth]{Figures/TempProfileClimatology_J5J7Z30.png}
    \caption{Globally and temporally averaged vertical temperature profile for for non-adaptive $2^\circ$ (J5), $1^\circ$ (J6) and $0.5^\circ$ horizontal resolution simulations compared with the adaptive simulation. (a) Full atmosphere. (b) Zoom near surface to show convergence with horizontal resolution.}
    \label{fig:TempProfileClimatologyJ5J7}
\end{figure}

\begin{figure*}[t]
\includegraphics[width=0.7\textwidth]{Figures/J5J7Z30_Temp_Vort_LonLat.png}
\caption{Longitude-latitude projections of temperature and vorticity in the surface layer from the adaptive {\sc wavetrisk}--simple physics simulation with seasons, with resolutions $2^\circ, 1^\circ, 0.5^\circ$ resolution, 10 soil layers and 30 atmosphere layer. Top: 5-year mean temperature. Middle: temperature on March~22 of year 5. Bottom: vorticity on March~22 of year 5, showing turbulent structure . Comparing this figure with the lower $2^\circ$ non-adaptive results~\ref{ClimatologyTempNonAdapt} shows how increasing the resolution leads to more intense vorticity filaments, but the temperature is relatively unchanged.}
\label{fig:ClimatologyTempVortAdapt}
\end{figure*}

The preceding statistical and instantaneousness results  confirm that the dynamically adaptive \ac{GCM} is indeed stable and able to capture the emergence of small scale turbulence and seasonal dynamics. Furthermore, the similarity of the high resolution adaptive climatology to the low resolution non-adaptive climatology suggests that adapting only on the prognostic variables is sufficient, at least for dry physics. The implicit effect of the physics on the dynamics is sufficient to resolve time- and scale-dependent dry physics phenomena. To be sure, however, we also compare the adaptive simulation results to an equivalent non-adaptive $0.5^\circ$ resolution simulation.



To get a more detailed idea of how well the adaptive simulation captures the weather dynamics, Figure~\ref{fig:Adapt_nonAdapt_Vorticity} compares instantaneous longitude-latitude projections of vorticity at height $P/P_S=0.5$ on 25 June of year 4, when the the weather is very different in the northern and southern hemispheres. The results are qualitatively similar, although the small scale structures in the northern hemisphere are not quite as well resolved in the adaptive simulation.  Note that we have set the relative tolerance to 2.5\%, and the tolerance could be decreased (or increased) to balance computational cost and fidelity of the weather as desired.  It is also important to note that because the flow is turbulent we cannot expect the details of the flow to match (e.g., shape and position of particular vortex filaments, number of eddies etc.). These figures also show that the turbulence during the summer solstice is extremely different in the two hemispheres, with a narrow band of band of small scale turbulence cenred around $30^\circ$~N and a large scale less intense turbulence extending over the entire southern hemisphere.  The adaptive grid shown in Figure~\ref{fig:Adaptive_grid} (bottom) shows that this difference in turbulence is well-captured by the adaptive dynamical grid which has essentially no refinement north of about $50^\circ$.
%
\begin{figure*}[t]
\includegraphics[width=0.7\textwidth]{Figures/Adapt_nonAdapt_Vorticity.png}
\caption{Vorticity on 25 June of year 4 at height $P/P_S=0.5$ (approximately layer 12). Top: adaptive with resolutions $2^\circ, 1^\circ, 0.5^\circ$. Bottom: non-adaptive with resolution $0.5^\circ$.}
\label{fig:Adapt_nonAdapt_Vorticity}
\end{figure*}
%

Taken together, the above results suggest that adapting on the prognostic dynamical variable alone is sufficiently accurate for the dry simple physics, at least up to $0.5^\circ$ resolution.  

A key feature of any \ac{GCM} is the energy budget, and in order to properly represent turbulent atmosphere dynamics we need to ensure that the energy spectrum of the adaptive simulation is similar to that of the finest non-adaptive simulation is important over all length scales. Because the wavelet multiscale approach adapts in both position and length scale, one of its strengths is its ability to capture the full turbulent energy spectrum with a much smaller number of computational elements than an equivalent spectral method~\citep{goldstein2005cvs_scales}.  Our goal now is to check that the adaptive simulation accurately captures the complete turbulent energy spectrum of the non-adaptive simulation. Since the flow is compressible, the rotational and divergent parts of the flow are both relevant. Recall that turbulence is characterized by a power law energy spectrum, where the power law $p, E(k)\propto k^p$ varies depending on the type of turbulence.

Figure \ref{fig:Spectra_scaling} shows the power laws for each vertical layer for the $0.5^\circ$ non-adaptive simulation. The slope $p$ is found from linear fits on a $\log-\log$ of energy spectrum for length scales between 600\,\unit{km} and 2000\,\unit{km}. The energy spectra were averaged over five years.  The power laws for the $0.5^\circ$ resolution simulations are between $-1.8$ and $-3$ for the divergent mode and between $-2.8$ and $-5.2$. For comparison, incompressible homogeneous isotropic three-dimensional turbulence has a power law of $-5/3$, while two-dimensional compressible turbulence has various spectra between $-2$ and $-4$. Overall, the divergent component has a shallower slope, meaning that it is more turbulent. Furthermore, the flow is the most turbulent in the atmospheric boundary layer for $P(z) > 600$\,\unit{hPa}. 

Figure~\ref{fig:Spectra} compares the rotational and divergent components of the energy spectra of the non-adaptive (left) and adaptive (right) simulation. Results are shown for layers 1 (993\,\unit{hPa}), 9 (713\,\unit{hPa}) and 27 (0.0111\,\unit{hPa}). The energy spectra for the adaptive simulation are similar to the non-adaptive spectra, but smoother. Most importantly the adaptive simulation captures the full range of active length scales and the correct turbulent slopes. The adaptive spectrum is smoother, especially at the smaller scales, due to the interpolation of the results to the finest grid resolution which is needed to compute the spherical harmonics energy spectra. Overall, even with the smoothening, the adaptive simulation captures the spectrum at all length scales even though it contains fewer degrees of freedom. 

\begin{figure*}[t]
\centering
\includegraphics[width=0.6\textwidth]{Figures/EnergyPowerLaws_J7Z30.png}
\caption{Energy spectra data averaged over year 5 for the non-adaptive $0.5^\circ$ resolution case: power laws for each layer obtained from linear fits for scales between 600\,\unit{km} and 2000\,\unit{km}. Note that the divergent spectra are much shallower than the rotational spectra and the flow is most turbulent at about 700\,\unit{hPa}.}
\label{fig:Spectra_scaling}
\end{figure*}

\begin{figure*}[t]
\centering
\includegraphics[width=0.49\textwidth]{Figures/EnergySpectraJ7Z30.png}
\includegraphics[width=0.49\textwidth]{Figures/EnergySpectra_J5J7Z30.png}
\caption{Energy spectra data averaged over year 5 for the $0.5^\circ$ resolution case. Rotational and divergent parts at layers 1 (993\,\unit{hPa}), 9 (713\,\unit{hPa}) and 27 (0.0111\,\unit{hPa}). Left non-adaptive simulation. Right: adaptive simulation with resolutions $0.5^\circ$, $1^\circ$, $2^\circ$. The adaptive simulation results have been interpolated to the finest resolution, $0.5^\circ$, for the spectral analysis.  The adaptive and non-adaptive results are consistent, although the adaptive results are smoother at small scales due to the interpolation.}
\label{fig:Spectra}
\end{figure*}

Finally, Figure~\ref{fig:Vorticity3D_J5J7Z30} shows isosurfaces of vorticity for the adaptive simulation on March 22 of year 5. The vorticity is primarily columnar, although there is some evidence of the intense atmospheric boundary layer structures visible in the eddy kinetic energy, eddy heat flux and eddy momentum flux shown in Figures~\ref{ZonalStatsAdapt_SpringFall}, \ref{ZonalStatsAdapt_SummerWinter}. The primarily vertical vortices suggests that adapting only horizontally should be fairly efficient for these climate dynamics.
%
\begin{figure*}[t]
\centering
\includegraphics[width=0.8\textwidth]{Figures/Vorticity3D_J5J7Z30.png}
\caption{Isosurfaces of vorticity on March 22 of year 5 ($\omega = \pm 10^{-4}$\,\unit{s^{-1}}) showing vertical structure of vorticity}
\label{fig:Vorticity3D_J5J7Z30}
\end{figure*}

Finally, we investigate whether a low resolution simulation can be successfully restarted from a saved checkpoint at much higher dynamically adaptive resolution.  This capability would allow investigation of climate dynamics at extremely high resolutions for shorter periods of time and quickly spinning up the model at low resolutions.  We restarted from a non-adaptive $2^\circ$ checkpoint at 100~days (June 30), with adaptive resolution increased by 8 times to $0.25^\circ$. This would correspond to a 64 times increase in horizontal resolution on a non-adaptive grid. We run the simulation for an additional 95 days until October 3.

Figure~\ref{fig:J5J8_compression} (top) shows the grid compression ratio as a function of time and figure~\ref{fig:J5J8_compression} (bottom) shows the adapted horizontal grid after 100 days. The restarted high resolution simulation ran stably with a very high average  grid compression ratio of 16.4. In comparison the adaptive simulation with maximum resolution $0.5^\circ$ had an average compression ratio of 2.3.  This result shows the potential of adaptivity to both achieve very him grid compression (and corresponding speed-up) and the ability to investigate very high horizontal resolutions for shorter periods. 
%
\begin{figure*}[t]
\centering
\includegraphics[width=0.6\textwidth]{Figures/J5J8Z30_compression.pdf}
\includegraphics[width=\textwidth]{Figures/J5J8Z30_adapted_grid.png}
\caption{Top: compression ratio for simulation restarted from a non-adaptive $2^\circ$ resolution checkpoint at 100~days (June 30) with maximum resolution $0.25^\circ$ (active scales are $2^\circ, 1^\circ, 0.5^\circ, 0.25^\circ$).  The average compression ratio is 16.4. Bottom: adapted grid after 100 days (October 3).}
\label{fig:J5J8_compression}
\end{figure*}


\conclusions\label{Conclusions}  
This paper presents the results of coupling a simple dry physics sub-grid scale model to the {\sc wavetrisk-atmosphere} dynamical core. We have four main goals:
\begin{enumerate}
\item 
Fully document the approximations and numerical implementation of the simple dry physics model introduced by \cite{Frederic_these}. 
\item 
Couple the simple physics package to the {\sc wavetrisk-atmosphere} multiscale adaptive dynamical core.  To do this, the simple physics code has been rewritten in FORTRAN~2008 and simplified.
\item 
Describe the 5-year climatology of the {\sc wavetrisk}--simple physics climate model, run non-adaptively at coarse resolution $2^\circ$ with 30 vertical atmosphere layers and 10 soil layers. The results without and with seasons are compared to an equivalent~\cite{HeldSuarez} climate model.
\item 
Explore an adaptive {\sc wavetrisk}--simple physics climate model coupling and assess the need for scale aware parameterizations or grid adaptation criteria. The coarsest resolution is $2^\circ$ and we allow two layers of refinement: $1^\circ$ and $0.5^\circ$, for a maximum grid compression ratio of 16 times. The adaptive results are compared with an equivalent non-adaptive $0.5^\circ$ resolution simulation.
\end{enumerate}
%
The basic approximations of the physical model are: dry atmosphere (i.e., no moisture), hydrostatic balance, compressibility and a floating vertical coordinate. The physics model is coupled to the dynamics as an implicit split step in time. It includes sub-grid scale parameterizations for radiation, vertical turbulent diffusion of velocity and potential temperature, a planetary boundary layer, surface fluxes, dry convection, diurnal cycle, seasons and, significantly, a soil column. The lack of moisture reduces model realism, but maintains simplicity and simplifies the interpretation of results. Model parameters, such as gravity and astronomical parameters, are set to Earth values.

Compared to the Held \& Suarez and the International Standard Atmosphere, the temperature profile for the simple physics--{\sc wavetrisk} climate model decreases is significantly warmer near the surface of the atmosphere and decreases monotonically with height. These differences are likely due to two limitations of the simple dry physics package, specifically the lack of moisture and the simplicity of the radiation parameterization (e.g., no ultraviolet absorption band in the stratosphere). We presented zonal statistics of temperature, zonal velocity, eddy kinetic energy, eddy heat flux, temperature variance and eddy momentum flux.  Without seasons, these statistics are qualitatively similar to those of the Held \& Suarez model, although the location of zonal jets is shifted and the eddy kinetic energy is about half as strong.  This suggests the additional slow time scales introduced by the soil column do not significantly affect the climatology. Unsurprisingly, adding seasons dramatically modifies all statistical quantities.

The adaptive simple physics--{\sc wavetrisk} climate model with horizontal resolutions $2^\circ, 1^\circ, 0.5^\circ$ is stable and accurate.  With relative tolerance 2.5\% produces grid compression ratios of between about 1.8 times (at the equinoxes) and 3.0 times (at the solstices).  This demonstrates the efficiency advantages of horizontal grid adaptivity, as it is able to track the active and inactive parts of the globe throughout the year. The zonal statistics are very similar to those of the $2\circ$ non-adaptive simulation and also agree well with the non-adaptive $0.5^\circ$ results.  

In addition to comparing the 5-year statistics, we also compared longitude-latitude projections of vorticity for the adaptive and non-adaptivity simulations on 25 June of year 4 when their is a strong north-south asymmetry in the turbulence.  The weather is characterized by a narrow band of small scale turbulence centered around $30^\circ$\,N and a large scale less intense turbulence extending over the entire southern hemisphere. The adaptive simulation successfully captures the large and small scale vortical structures in each hemisphere.  Finally, we presented isosurfaces of vorticity for March 22 of year 5. The vorticity is primarily columnar, but also includes intense atmospheric boundary layer structures, which are also visible in the statistics.

To complement the physical space comparisons, we compared energy spectra for the rotational and divergent components for adaptive and non-adaptive simulations. These results confirm that the adaptive simulation accurately captures the full range of active turbulent scales, although there is some smoothing due to interpolation of the adaptive grid onto the full grid (needed to compute the spherical harmonics transform).  Vertical profiles of the energy spectrum power law show strong dependence on depth, with the most turbulent region in the boundary layer, $P(z)>600$\,\unit{hPa}.

Overall, we conclude that the simple dry physics model can be stably coupled to an adaptive dynamical core.  Scale-aware parameterizations are not necessary, at least for dry physics. 

The climatology presented here could provide the basis for new dynamical core test cases since the simple dry physics can be easily coupled to different dynamical cores.  The simple physics includes relatively sophisticated radiation, turbulence, convection and diurnal/seasonal forcing and thus fills a gap between the Held \& Suarez model and realistic moist physics models. Future work could explore an intercomparison of a set of dynamical cores coupled to the same simple dry physics package. 

Although scale-aware physics parameterization is not necessary for this dry physics model, it is quite likely to be necessary for moist physics.  In future work, we will investigate coupling the adaptive {\sc wavetrisk} model to a basic moist physics, such as~\cite{kessler1969water}. This model has already been coupled to the Held \& Suarez model~\cite{ullrich2016moistheld}. 

We hope that the results presented here will stimulate interest in~\cite{Frederic_these}'s interesting dry physics model and has demonstrated the potential of adaptive dynamical cores to improve the efficiency and accuracy of climate modeling.  Our results may also help climate modelers simplify the tuning required when grid resolution is increased by focusing their attention on sub-models not included in the simple dry physics.

\codedataavailability{The {\sc wavetrisk} code, including the simple dry physics and Held--Suarez physics, are available at 
\url{https://github.com/kevlahan/wavetrisk_hydrostatic}. 
The test case {\tt ~/test/climate} is used to run the simulations. The data presented in this paper are available at this globus link:
\href{https://globus.alliancecan.ca/file-manager?origin_id=2cd3329f-064f-4489-b634-eeeafa59f1a7\&origin_path=\%2F}{Globus}
.} 


\appendix
\section{Software Engineering Design Decisions and Challenges\label{sec:software}}    %% Appendix A  
To couple the simple dry physics package to the dynamical core {\sc wavetrisk} certain design decisions were made in order to facilitate the physics-dynamics coupling and overcome challenges encountered due to dynamic adaptivity. This appendix section will summarize the general steps needed to couple the simple dry physics package to a dynamical core and the decisions made to couple the package to {\sc wavetrisk}.

In order to utilize the physics package as a test case and couple it with a dynamical core an interface is used. This interface is the connector and facilitator of the entire model and will vary among climate models due to different dynamical core architectures and their methods of incorporating test cases. In terms of coupling the physics to {\sc wavetrisk}, we made the design decision to create an interface formatted as one of {\sc wavetrisk}'s test cases. The interface calls all physics initialization routines through a wrapper function which includes setting up the plugins of the physics and reading all constants and constraints set by the dynamics. As the facilitator, the interface couples the dynamics and physics by calling their associated timestep routines, usually in a split step. When calling the physics timestep a few aspects need to be taken into account, the precision of the both the dynamics and physics, the format of the main prognostic variables stored and the data structures of the prognostic variables, which is highly dependent on the grid. 

The physics package is written in single precision, using only 9 digits to represent floating point real values, while {\sc wavetrisk} is written in double precision. To avoid converting to and from both data types, before and after the physics step, we decided to compile the physics package with the Fortran flag, {\tt -freal-4-real-8}. This flag converts all real numbers (single precision) in the physics package to real (8) values (double precision), allowing for compatibility with {\sc wavetrisk}.

The key input variables required for a physics package time step includes temperature, zonal velocity and meridional velocity at each vertical layer, along with the pressure and geopotential of the vertical columns. Calculating the pressure and geopotential requires the simple integration from the surface upwards for each column and the utilization of the hydrostatic approximation. In comparison, the temperature and velocities require a more in-depth conversion. {\sc wavetrisk} saves the mass weighted potential temperature as its prognostic variable. Therefore, when calling the physics, the interface converts the value, at each vertical layer of a column, to potential temperature and then to the required temperature input, in Kelvin, of the physics. Upon return from the physics call, the interface converts the temperature tendencies back to mass weighted potential temperature using the reverse process. Regarding the zonal and meridional velocities, wavetrisk does not store the velocities at the center of the node, instead it stores three velocity edge values. To retrieve the zonal and meridional velocities at the node, interpolation using the three edge values and neighbouring columns' edge values is required. The inverse procedure is also applied when the physics returns the zonal and meridional tendencies for each column. 

In general, the data structures storing the prognostic variables of both the dynamics and physics are different. {\sc wavetrisk} utilizes a hybrid data structure which is optimal for the grid chosen and for adaptivity and load balancing associated with the use of parallelism (see~\citep{wavetrisk_datastruct} for more on the {\sc wavetrisk} grid architecture and data structure). Conversely, the physics package utilizes a regular data structure to store the values for each column. The regular data structure is a 2D array, with each row representing a different vertical column and each column of the array represents a vertical layer. Upon initialization of the physics package, the number of vertical columns to be sent at each time step is fixed. At each physics step, the package takes in the fixed number of vertical columns, the time step, the simulation day, the fraction of the day and the columns' prognostic variable values, pressure and geopotential at each vertical layer. A major overhead of the simulation, even in the case of parallelism, is the retrieval and conversion of all vertical columns from the hybrid data structure to the required regular structure for the input of the prognostic variables, as well as the conversion of the output tendencies from the regular to hybrid data structure. 

The use of a 2D data structure constrains the package as it expects the user to call the physics for all vertical columns at once. However, another major obstacle this structure produces is in the case of parallelism and adaptivity. {\sc wavetrisk} utilizes {\tt mpi} and domain decomposition of the columns of the sphere. When adaptivity is used, a major concern is ensuring correct load balancing on each CPU. Therefore, {\sc wavetrisk} allows for the re-balancing of the domains on each CPU, to ensure approximately equal computational loads. The complication arising with the ability of rebalancing is that the physics package saves the soil and surface temperatures internally for each column. To reduce the overhead and overcome the rebalancing obstacle, we have modified the physics step to call the physics for vertical columns individually and the step is completed when the last column's tendency returns from the call. However, to allow for single column transmission, the physics package and {\sc wavetrisk} required some fine tuning and additions.

The major addition to the physics package is the creation of the single column module. It allows the model to call the physics package for each column individually by providing the soil and surface temperatures to the dynamics. This module contains a wrapper subroutine that the interface calls, including all required inputs for the original physics call plus the added soil and surface temperatures. Once the wrapper subroutine is called, it updates the soil and surface temperatures for the column and calls the original physics routine. The subroutine outputs the original tendencies and the new updated soil and surface temperatures for the interface to save. The module also contains a subroutine that allows for the change in latitude and longitude required when sending a different column at each physics call. A major change to the physics package is the input of the number of desired soil layers when initializing the parameters of the grid in the package. Originally this was fixed to 10, but a variable number of layers is desired to test the importance of the soil model. 

The only major change to {\sc wavetrisk} was in the grid initialization. To send single columns to the physics with the possibility of rebalancing, {\sc wavetrisk} must save the soil temperatures and surface temperature of each column. The vertical layers were therefore extended, below the surface, employing  negative and zero indexes for each vertical column to represent the soil layers ($k<0$) and surface ($k=0$) temperature respectively. The interface facilitates this by setting the number of soil layers, before all grids are initialized. The only drawback of extending the grid is that all prognostic variables' data structures needs to be extended, even though only the mass-weighted potential temperature negative indexes are used to save the soil and surface values. Therefore, there is additional unnecessary memory overhead. In future, it is recommended that only the required variable's data structure be extended, but this will require a more extensive update, not only to the architecture of the grid, but also to the architecture of the calls for the dynamics step. 

Another {\sc wavetrisk} design decision was the format and storage of the extra temperature values. It was decided that the soil and surface temperatures would be stored in the mass-weighted potential temperature's data structure in the negative and 0 indexes respectively. More importantly, the temperatures would not be converted to mass-weighted potential temperature, but would be retained as temperature (units of Kelvin). This decision was made to avoid the extra overhead that comes with the conversion, as the values are not used by the dynamics.

Finally, while the main decision to send individual grid columns to the physics required extra temperature values to be stored by the dynamics, no other physics variables need to be stored. In many models, as in the {\sc wavetrisk-ocean} model configuration, a new dynamics variable tke is required to keep track of the turbulent kinetic and the tke evolves over time with the physics' k-epsilon model. However, since the simple physics model's vertical diffusion is a first order closure model, no further variables need to be stored and the interpolation of any internal physics variables onto the adaptive grid is not required.

\subsection{The workflow of the interface and physics call}
%
To illustrate the sequence of steps taken by the interface Fig.~\ref{InterfaceProgram} displays the workflow of the interface written to couple {\sc wavetrisk} to the physics package. Furthermore, to illustrate the physics-dynamics coupling challenges mentioned above, Fig.~\ref{PhysicsStepWorkflow} displays a flowchart of the algorithm used by the interface and required for the a single physics (split) step. The steps in black are the extra steps required to change prognostic variables to the required form, for example mass-weighted potential temperature to temperature or vice versa. The yellow step represents the call to the wrapper subroutine in the single column module added to physics package to allow the architecture of sending individual columns. 

\begin{figure}[ht]
\centering
\includegraphics[scale=.27]{Figures/Interface_Program.png}
\caption[Interface Program Work Flow]{A flow chart of the interface (physics-dynamics) coupling the physics and dynamics. The workflow shows the the steps of the simulation until the end.}
\label{InterfaceProgram}
\end{figure}

\begin{figure}[ht]
\centering
\includegraphics[scale=.2]{Figures/Physics_SplitStep_Workflow.png}
\caption[Workflow of the {\sc wavetrisk} Physics Split Step Call]{A flow chart of the physics step. Steps in black contain overhead concerning the conversion of the prognostic variables from one type to another (i.e.\ mass-weighted potential temperature to temperature ($K$)). The yellow process represents the physics call of the added single column module of the physics.}
\label{PhysicsStepWorkflow}
\end{figure}
\noappendix      

\authorcontribution{F Hourdin developed the original simple physics model in his PhD thesis. GCJ, NKRK prepared the paper with contributions from TD. NKRK developed the model code with contributions from GCJ and TD for the dry physics module. GCJ and NKRK performed the model simulations using {\sc wavetrisk-atmosphere}.}

\competinginterests{The contact author has declared that none of the authors has any competing interests.}

\begin{acknowledgements}
NKRK and GCJ were supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) under grant no. RGPIN-2024-05282. 
\end{acknowledgements}

\bibliographystyle{copernicus}
\bibliography{bibliography}

\end{document}
