remarkRemark exampleExample propProposition corCorollary Well-posedness and
regularity at cloud edgeA. Remond-Tiedrez, L. Smith, and S. Stechmann

#  A nonlinear elliptic PDE from atmospheric science: well-posedness and
regularity at cloud edge

Antoine Remond-Tiedrez () Department of Mathematics, University of
Wisconsin–Madison, Madison, WI, USA and Department of Applied Mathematical and
Theoretical Physics, University of Cambridge, Cambridge, UK
<EMAIL_ADDRESS>Leslie M. Smith () Department of Mathematics and
Department of Engineering Physics, University of Wisconsin–Madison, Madison,
WI, USA<EMAIL_ADDRESS>Samuel N. Stechmann () Department of Mathematics
and Department of Atmospheric and Oceanic Sciences, University of
Wisconsin–Madison, Madison, WI, USA<EMAIL_ADDRESS>
###### Abstract

The precipitating quasi-geostrophic equations go beyond the (dry) quasi-
geostrophic equations by incorporating the effects of moisture. This means
that both precipitation and phase changes between a water-vapour phase
(outside a cloud) and a water-vapour-plus-liquid phase (inside a cloud) are
taken into account. In the dry case, provided that a Laplace equation is
inverted, the quasi-geostrophic equations may be formulated as a nonlocal
transport equation for a single scalar variable (the potential vorticity). In
the case of the precipitating quasi-geostrophic equations, inverting the
Laplacian is replaced by a more challenging adversary known as potential-
vorticity-and-moisture inversion. The PDE to invert is nonlinear and piecewise
elliptic with jumps in its coefficients across the cloud edge. However, its
global ellipticity is a priori unclear due to the dependence of the phase
boundary on the unknown itself. This is a free boundary problem where the
location of the cloud edge is one of the unknowns. Potential-vorticity-and-
moisture inversion differs from oft-studied free-boundary problems in the
literature, and here we present the first rigorous analysis of this PDE.
Introducing a variational formulation of the inversion, we use tools from the
calculus of variations to deduce existence and uniqueness of solutions and de
Giorgi techniques to obtain sharp regularity results. In particular, we show
that the gradient of the unknown pressure or streamfunction is Hölder
continuous at cloud edge.

###### keywords:

potential vorticity inversion, quasi-geostrophic equations, moist atmosphere,
free-boundary problems

Primary: 49J10, 35B65, 86A10; Secondary: 35J20, 35Q86, 76U60.

Table of contents

  1. 1.

Introduction 1

  2. 2.

Discussion 2

  3. 3.

Background 3

  4. 4.

Variational formulation and basic properties of the energy 4

  5. 5.

Well-posedness 5

  6. 6.

Smoothing 6

Appendix

  1. A.

Tools from convex analysis A

  2. B.

Tools from regularity theory B

  3. C.

Rosetta stone C

  4. D.

Comparison of the variational and conserved energies D

Note to the reader: Section 1 acts as a “shortest path” to the statements of
the main result of this paper. Section 2 discusses the main difficulties
encountered, provides a map of the remainder of the paper and its arguments,
and situates the results herein with respect to the broader literature. More
details of the model itself are presented in Section 3 since the model is not
well-known to the mathematical community. Finally Sections 4–6 are the main
body of the paper and contain the arguments culminating in the central
theorems stated in Section 1 below.

##  1 Introduction and statement of the results

An important success story in atmospheric dynamics is the description of
incompressible fluids in the limit of fast rotation and strong stratification.
This limit is relevant on the synoptic scale111 The word “synoptic” has Greek
roots: “syn–” means “together” and “–optic” means “seen”. The synoptic scale
is thus that at which a “general picture” emerges. In concrete terms it
usually refers to a spatial scale of about 1,000 km and a temporal scale of a
few days.  for the description of the troposphere and was introduced by
Charney in [Cha48] as the _quasi-geostrophic_ (QG) model [Maj03, Val17]. One
of its key successes is that it reduces the dynamics to a single scalar
variable, namely the _potential vorticity_ which we will denote as $PV$
[Ert42]. Crucially for our discussion here: phrasing the entire QG dynamics in
terms of the $PV$ requires the inversion of a Laplacian – this is known as
potential vorticity inversion, or _$PV$ inversion_ [HMR85]. This framework is
similar to the vorticity formulation of the incompressible Euler equations
where the Biot-Savart law requires the inversion of a Laplacian.

The classical QG model does not take into account moisture, despite the key
role that moisture plays in atmospheric dynamics, most notably because water
underpins important mechanisms for energy transfer in the atmosphere. As water
evaporates it absorbs heat, while as water condensates it releases the latent
heat it had stored. This motivates consideration of _moist_ variants of the QG
system.

In this paper we will consider one such moist model: the _precipitating quasi-
geostrophic_ (PQG) model introduced in [SS17]. In that model the dynamics are
now fully described by two variables: an _equivalent_ potential vorticity
$PV_{e}$ which differs slightly from the one discussed above (more on that in
Section 3) and a moist scalar variable $M$. By contrast with the classical
case in which moisture is not accounted for and where $PV$ inversion requires
the inversion of a Laplacian, when moisture is involved we now have to deal
with $PV_{e}$-and-$M$ inversion, which relies on the inversion of a more
complicated operator. The complications involve both nonlinearities and lack
of smoothness at cloud edge. The object of this paper is the first rigorous
analysis of this operator and of $PV_{e}$-and-$M$ inversion.

We need to discuss one more aspect of the model before we can write down
$PV_{e}$-and-$M$ inversion, namely the manner in which the moisture is
incorporated in the model. The moisture is included in the model via a single
scalar unknown $q$ which describes the total water content (water vapour plus
condensed water). For simplicity we will assume without loss of generality
that saturation occurs at $q=0$.

Since different physics will occur depending on the sign of $q$, we introduce
the following notation for convenience. The unsaturated phase
$\left\\{q<0\right\\}$ containing water vapour only is indicated by the
Heaviside function $H_{u}=\mathds{1}(q<0)$. For the saturated phase
$\left\\{q\geqslant 0\right\\}$, which contains water in the forms of both
vapour and liquid, we write $H_{s}=1-H_{u}=\mathds{1}(q\geqslant 0)$ (see also
Figure 1).

Figure 1:  A pictorial depiction of the interface separating the two phases of
water considered here (recall that, without loss of generality, saturation is
taken to occur at $q=0$). The blue region on the left corresponds to the
unsaturated region where $H_{u}=1$ (the outside of a cloud) whereas the pink
region on the left where $H_{s}=1$ corresponds to the saturated region (the
inside of a cloud). The level set $\left\\{q=0\right\\}$ is the cloud edge.
Recall that $H_{u}=\mathds{1}(q<0)$ and $H_{s}=\mathds{1}(q\leqslant 0)$ such
that $H_{u}+H_{s}=1$.

$PV_{e}$-and-$M$ inversion then takes the following form,

$\Delta_{h}p+\partial_{3}\left(\frac{1}{2}\left(M+\partial_{3}p\right)H_{u}+\left(\partial_{3}p\right)H_{s}\right)=PV_{e},$
(1)

where $PV_{e}$ and $M$ are given and we seek to solve for the pressure $p$.222
In the precipitating quasi-geostrophic system the pressure acts as a
streamfunction (as in the dry case). Solving for the pressure thus means
solving for the velocity. Since the velocity transports the active scalars
$PV_{e}$ and $M$, which in turn determine the pressure through
$PV_{e}$-and-$M$ inversion, the system is closed. This is discussed in more
detail in Section 3.  Note that throughout this paper we will work on the
$3$-torus $\mathbb{T}^{3},$ and we use the notation
$\Delta_{h}=\partial_{1}^{2}+\partial_{2}^{2}$, where $x_{3}$ denotes the
direction of both rotation and stratification (this is the so-called
“traditional approximation” [Eck60]). The PDE (1) is nonlinear since $H_{u}$
and $H_{s}$ are determined by the water content $q$, which is in turn a
function of $M$ and $\partial_{3}p$ (see Section 3 for a more detailed
discussion of the precipitating quasi-geostrophic model, which includes a
derivation of (1) and more details on how $M$ and $\partial_{3}p$ determine
$q$). Immediately, (1) introduces two key questions.

Question 1: Is equation (1) truly elliptic, and thus solvable? Certainly it is
elliptic in each of the saturated and unsaturated phases, but how does the
dependence of the phase boundary on the unknown $p$ impact the overall
ellipticity and solvability of the problem?

A first attempt at answering this would seek to smooth out the Heavisides
appearing in (1). This raises subtle issues discussed in Remark 6.1.

Question 2: What quantities, if any, are guaranteed to be continuous across
the interface between the two phases? Conversely, are there examples for which
certain quantities jump across the interface?

The fundamental observation required to answer these questions is the
following. We may reformulate $PV_{e}$-and-$M$ inversion by invoking the
minimum function instead of explicitly using Heavisides (this is discussed in
Section 3 below) to obtain

$\Delta
p+\frac{1}{2}\partial_{3}\left(\min\left(M-\partial_{3}p,\,0\right)\right)=PV_{e}.$
(2)

We will use the notation involving the minimum function throughout this paper
but note that we could alternatively have used the fact that
$\min(f,\,0)=f^{-}$, where $f^{-}$ denotes the negative part of a function
$f$, to phrase $PV_{e}$-and-$M$ inversion in terms of the negative part of
$M-\partial_{3}p$. The formulation (2) of $PV_{e}$-and-$M$ inversion allows us
to answer both of the key questions above, as shown by the well-posedness
result stated below. Note that all functional spaces are assumed to be taken
over the torus $\mathbb{T}^{3}$ unless specified otherwise. Moreover the main
spaces of interest will be
$L^{2}\vcentcolon=\left\\{u:\int_{\mathbb{T}^{3}}u^{2}<\infty\right\\}$,
$\mathring{H}^{1}~{}\vcentcolon=~{}\left\\{u\in L^{2}:\nabla u\in L^{2}\text{
and }u\text{ has average zero}\right\\}$, and its dual
$H^{-1}\vcentcolon={\big{(}\mathring{H}^{1}\big{)}}^{*}$.

######  Theorem 1.1.

For every $M\in L^{2}$ and $PV_{e}\in H^{-1}$ there exists a unique weak
solution $p\in\mathring{H}^{1}$ of (1). Moreover if $PV_{e}$ and $M$ are
smooth then

  1. 1.

$p\in C^{1,\alpha}$ globally, including across the interface, in the sense
that the gradient of $p$ is Hölder continuous with exponent $\alpha$ and

  2. 2.

both the unsaturated and saturated phase, $\left\\{q<0\right\\}$ and
$\left\\{q>0\right\\}$ respectively, are open sets and $p$ is smooth in each
phase.

Finally: this regularity result is sharp in the sense that there exist smooth
$PV_{e}$ and $M$ for which the unique solution $p$ of (1) belongs to
$C^{1,\,1}$ but not $C^{2}$, i.e. the gradient of $p$ is Lipschitz but not
differentiable.

##  2 Discussion: difficulties, overview, and related literature

The main difficulty in the study of $PV_{e}$-and-$M$ inversion in its native
form (1) is the presence of the Heavisides. From a modelling perspective it is
remarkably useful to invoke the Heavisides $H_{u}$ and $H_{s}$. They allow us
to clearly see how physical relations depend on water according to the phase
in which they occur. From a mathematical perspective, however, they introduce
an ouroboros-esque danger of eating your own tail. Given fixed Heavisides, (1)
is piecewise elliptic and so its solvability seems all but assured. However,
these Heavisides are not fixed and depend on the solution. Here lies the
difficulty: the solution is determined by the Heavisides, which themselves are
determined by the solution. The snake is eating its own tail. Note that this
is not a surprising difficulty to encounter since this is the typical issue
that arises when tackling free boundary problems.

On the other hand, rewriting (1) in the form of (2) allows for progress
because it folds the free boundary $\left\\{M=\partial_{3}p\right\\}$
explicitly into the nonlinearity. This means that solving the PDE requires no
_a priori_ knowledge of the free boundary: we may simply find the solution via
a variational characterization and _then_ determine the location of the free
boundary from the solution. The snake’s tail is out of its mouth.

To be more precise, reformulating $PV_{e}$-and-$M$ inversion in the form of
(2) helps us in two ways. First, the formulation (2) leads more directly to a
variational formulation of $PV_{e}$-and-$M$ inversion. The Direct Method of
the Calculus of Variations may then be brought to bear to deduce the first
part of Theorem 1.1, as shown in Section 5, culminating in Theorem 5.1. In
particular, we provide one way of answering our initial question asked
immediately after first encountering $PV_{e}$-and-$M$ inversion in its form
(1) involving Heavisides: demonstration that the energy is strongly convex
(see Proposition 1) verifies that indeed $PV_{e}$-and-$M$ inversion is an
elliptic problem.

Note that there is another candidate for the energy characterising
$PV_{e}$-and-$M$ inversion. This alternate candidate is arguably more natural
since it corresponds to the energy conserved by the PQG dynamics. However, for
reasons detailed in Remark D.4, this conserved energy is not suitable and thus
highlights the care required to identify an appropriate variational principle.

The second way in which the formulation (2) helps us is with regard to
regularity. Indeed, while differentiating (1) would prove challenging due to
the presence of Heavisides, differentiating that same equation in the form of
(2) is now straightforward. It tells us that $u=\partial_{i}p$ solves

$\partial_{i}PV_{e}=\Delta
u-\frac{1}{2}\partial_{3}\left(\mathds{1}(M<\partial_{3}p)\partial_{3}u\right),$
(3)

where we have used the notation ${(v\otimes w)}_{ij}=v_{i}w_{j}$ for any two
vectors $v$ and $w$. In particular, since the coefficient matrix
$A=I-\frac{1}{2}H_{u}e_{3}\otimes e_{3}$ is uniformly elliptic, we may use the
tools from de Giorgi’s regularity theory to deduce the sharp regularity result
for $PV_{e}$-and-$M$ inversion recorded in the second part of Theorem 1.1. We
note that the regularity theory presented in Section 5 is classical albeit
difficult to find in detail in the literature. References [Vas16] and [CV10]
are excellent introductions to the simplest setting in which de Giorgi tools
bear fruition, when there is no forcing, and [HL11] contains a version of the
regularity theory presented here but without proof. The argument presented in
detail here emphasizes how to adapt the classical de Giorgi result to the case
of nonzero forcing through the use of the Campanato characterisation of Hölder
continuity. We hope that this presentation can serve the community by filling
that expository gap in the literature. Note also that (3) provides another way
of answering our first question of whether or not $PV_{e}$-and-$M$ inversion
is truly elliptic: the ellipticity of the coefficient matrix $A$ discussed
above answers that question in the affirmative.

In the last section of this paper we provide some computations on how to
smooth out $PV_{e}$-and-$M$ inversion. These computations help understand why
analysing the ellipticity of $PV_{e}$-and-$M$ inversion directly from (1) is
quite subtle, and they are also of independent interest. Indeed, they provide
fertile ground for further analysis: when solving $PV_{e}$-and-$M$ inversion
numerically, it may be beneficial to first solve the smoothed-out version, and
then let the smoothing parameter decrease to zero. It also remains to be seen
whether the smoothed out interfaces are close to the true interface – this
would be guaranteed if we had $C^{1,\,\alpha}$ estimates on the distance
between the smoothed out minimizers and the true minimizer. As a first step we
provide $H^{1}$ estimates here.

Ultimately, this paper is also intended to bring the attention of the PDE
community to a free boundary problem not yet studied. At the end of this
section we discuss how this problem differs from those most commonly studied
in the literature. The well-posedness result provided here is but the first
step in the more careful analysis of this particular free boundary problem,
and also serves to give confidence in numerical studies of the PQG system that
rely on $PV_{e}$-and-$M$ inversion (such as [HESS21]). Moreover, important
open questions remain. For instance we do not yet know anything about the
regularity of the interface, even conditionally based on assumptions on the
regularity and/or compatibility of $PV_{e}$ and $M$.

The rigorous analysis of moist atmospheric models has only recently attracted
the attention of the PDE community. As we review some of that work, we focus
on the rigorous treatment of models involving phase transitions (some earlier
works studied moist quantities that were purely advected by the velocity
field, without phase changes).

The first rigorous foray into moist atmospheric dynamics with phase changes
can be traced back to a series of work [CZT12, CZFTT13, BCZT14] which studies
a two-phase (water vapour and condensed water) model where the velocity field
is taken to be given. The velocity is then evolved dynamically via the
primitive equations in [CZHK+15] and topography is considered in [LM20]. The
same model, but with non-constant water saturation ratio, is studied in [TW15,
TW16].

There are then two branches of study for _three_ –phase models (which involve
water vapour and distinguish two kinds of condensed water: cloud water and
precipitating water).333Note that the version of the PQG model used here is a
two-phase version in which the condensed water phase can be cloud water that
does not fall, or precipitating water (rain) with prescribed fall speed, but
the version here does not consider co-existence of clouds and rain or the
conversion from cloud water to rain water. Other versions of precipitating
quasi-geostrophic models, with three phases of water or other additional
complexities in cloud microphysics, can also be defined; see discussion in
[SS17, WSSM19, WSS+20]. One branch relies on cloud microphysics from [Gra98],
beginning with [CHT+18] where the velocity is given. In [TL22] the velocity
then solves the primitive equations. Another branch relies on cloud
microphysics from [MK06], beginning with [HKLT17] where the velocity is given.
In [HKLT20] the velocity solves the primitive equations instead.

The incorporation of additional phases due to ice is discussed in [CJTT21]. By
contrast, earlier works which omit cloud ice are thus sometimes referred to as
“warm cloud” models.

Note that all works mentioned above either treat the velocity as prescribed or
as governed by the primitive equations. Following [CT07], a well-understood
framework relying on the barotropic–baroclinic decomposition of the velocity
is in place for the analysis of the primitve equations and its extensions. To
the best of the authors’ knowledge, there are as of yet no rigorous
mathematical treatments of moist models leveraging other dynamical models for
the velocity field. One such dynamical model is the precipitating quasi-
geostrophic system discussed here. In order for this system to be well-
formulated as a nonlocal transport system, $PV_{e}$-and-$M$ inversion needs to
be understood first. We thus provide the initial steps toward rigorous
analysis of a moist dynamical core related to, but different from, the
primitive equations.

Free boundary problems have by contrast garnered significant interest from the
mathematical community. Rather than provide an extensive review of the
literature, we direct the reader to [CSV15, FS15] and references therein.
Instead, we only highlight the key difference (and thus sources of interest)
of $PV_{e}$-and-$M$ inversion compared with oft-studied free boundary problems
in the literature. Here the free boundary depends on the _gradient_ of the
unknown, and not on the unknown itself. The latter case is the more typical
free boundary problem known as the obstacle problem (see the book [PSU12]).

Note that the dependence of the free boundary on the gradient is not unseen in
the literature – it is studied for example as “Model Problem B”, one of three
instances of standard classes of free boundary problems studied in [PSU12].
That being said, this “Model Problem B” is isotropic since the free boundary
is the boundary of the set $\left\\{\lvert\nabla u\rvert>0\right\\}$, whereas
the problem of interest here is _anisotropic_ since the _vertical_ derivative
of $p$ determines the free boundary. It is also worth noting that, although
not unseen in the literature, gradient-dependent problems are less prominent.

##  3 Background and derivation of $PV_{e}$-and-$M$ inversion

In this section we discuss the PQG equations in more detail since that model
is recent and not yet well-known in the mathematical community. In particular
we provide a derivation of $PV_{e}$-and-$M$ inversion in both its form (1) and
(2).

In order to motivate the form taken by $PV_{e}$-and-$M$ inversion, we will
first discuss how $PV$ inversion arises in the dry case. We will then see how
these arguments are modified once moisture is incorporated as per the PQG
model introduced in [SS17]. This will enable us to make the form of
$PV_{e}$-and-$M$ inversion precise and to state the main results we obtain
pertaining to the solvability of that inversion.

The dry case. We now turn our attention to a discussion of the dry QG model.
While the complete dynamics may be characterised solely in terms of the
potential vorticity $PV,$ it is helpful to describe the system by introducing
the auxiliary unknowns: the velocity $u$, the potential temperature444 The
potential temperature is a quantity defined in terms of the temperature and
the (compressible) pressure. It is often more convenient to work with
potential temperature than temperature itself (as is the case for dry QG)
since, unlike the temperature, the _potential_ temperature is purely advected
by the flow.  $\theta$, and the pressure $p$.

With these auxiliary variables in hand we can view the quasi-geostrophic
system as a combination of three features. The first feature is that the
potential vorticity defined as $PV={(\nabla\times u)}_{3}+\partial_{3}\theta$
is purely advected, i.e. $\left(\partial_{t}+u\cdot\nabla\right)PV=0$. The
second feature is that the velocity is determined by _geostrophic balance_
$u_{h}=\nabla_{h}^{\perp}p,$ which arises from the effect of the Coriolis
force and assuming fast rotation. Here we have introduced the notation
$v_{h}=(v_{1},\,v_{2},0)$ to denote the horizontal components of any 3-vector
$v,$ as well as the notation $w^{\perp}=(-w_{2},\,w_{1})$ to denote the
$\frac{\pi}{2}$–rotation of any 2-vector $w$. Then we write
$v_{h}^{\perp}=(-v_{2},\,v_{1},\,0)$. The third feature is that the potential
temperature is determined by _hydrostatic balance_ $\theta=\partial_{3}p,$
which arises due to strong stratification. Plugging these two balances into
the definition of the potential vorticity tells us that $PV=\Delta p,$ and
solving this elliptic PDE is known as _$PV$ inversion_.

Noting that the velocity transporting the $PV$ may be reconstructed from the
$PV$ itself via $u=\nabla_{h}^{\perp}\Delta^{-1}PV,$ one can see that the QG
system is a nonlinear and nonlocal transport equation for the potential
vorticity. This is similar to how the vorticity form of the incompressible
Euler equations is a nonlinear and nonlocal transport equation for the
vorticity.

The moist case. The question is now: “What happens when moisture is
incorporated into the model?” In order to answer that question, we first
recall that the moisture is modeled via a single scalar unknown $q$ which
describes the total water mixing ratio, and that for simplicity the saturation
is taken to occur at $q=0$. To be precise: $q$ is the mixing ratio measured as
kilograms of total water per kilogram of dry air [MK06, HDMSS13]. We then
introduced $H_{u}=\mathds{1}(q<0)$ and $H_{s}=1-H_{u}=\mathds{1}(q\geqslant
0)$ to make it easier to see how physical relationships vary depending on
whether they occur in the unsaturated phase $\left\\{q<0\right\\}$ or in the
saturated phase $\left\\{q>0\right\\}$.

While geostrophic balance remains the same once moisture is considered,
hydrostatic balance only appears unaffected at first glance. Indeed, it is
still true that $\partial_{3}p=\theta$, but it is now essential to split the
potential temperature $\theta$ into two pieces, writing
$\theta=\theta_{e}-qH_{u}$. Here $\theta_{e}$ denotes the _equivalent_
potential temperature, and corresponds to the potential temperature that a
parcel of air would have if it were heated by converting its water vapour
content into condensed water, and $qH_{u}$ accounts for the latent heat stored
in water vapour. Once moisture is incorporated into the model, $\theta$ is no
longer purely advected as was the case for dry dynamics, since source terms
arise due to condensation and evaporation. It is therefore $\theta_{e}$ which
is purely advected, and thus the variable used to define moist potential
vorticity $PV_{e}$. Finally, the moist potential vorticity is defined as
$PV_{e}={(\nabla\times u)}_{3}+\partial_{3}\theta_{e}$ (using $\theta_{e}$ and
not $\theta$ as was done in the dry case).

Physically speaking, it is important to introduce the concept of a slowly
varying quantity whose dynamics is appropriate to describe temporal and
spatial variations on synoptic scales. In the atmospheric science literature,
slowly varying quantities are also called “balanced” or simply “slow.”
Mathematically, identification of appropriate slow variables follows from a
distinguished limiting process of the governing equations, in this case either
the dry or moist rotating Boussinesq equations. While dry $PV$ is slowly
varying in the context of dry dynamics, it is moist $PV_{e}$ that is slowly
varying in the presence of moisture and phase changes [SS17]. By contrast, the
variables $u$, $\theta_{e}$ and $q$ exhibit wave-like behavior, generally
associated with smaller spatial scales and faster temporal scales. As slow
variables in their respective dry and moist environments, $PV$ and $PV_{e}$
are of primary importance to meteorology because their dynamics is associated
with synoptic-scale weather patterns.

The description of moist synoptic-scale dynamics requires at least two slowly
varying unknowns, instead of of the single $PV$ variable in the dry case. This
is because there is a new unknown quantity $q$, without any new balances
appearing (which would reduce the number of unknowns needed to fully describe
the system). Following from a distinguished limiting analysis of the two-phase
moist Boussinesq equations, the new slow variable is denoted $M$ and defined
as $M=\theta_{e}+q$ [SS17].

Note that we have, in the previous discussion and throughout the remainder of
this paper, set a slew of physical constants to unity. A “Rosetta stone” of
sorts is provided in Appendix C, which discusses the PQG model when these
constants are present in their fully glory. Appendix C acts as a dictionary
between the notation of [SS17] and this paper, and makes clear the fact that
that there is no loss of generality in considering the simpler looking
$PV_{e}$-and-$M$ inversion (1).

Derivation of $PV_{e}$-and-$M$ inversion. We are now ready to derive (1) and
(2). Using hydrostatic balance and the moist variable $M$ allows us to write
the equivalent potential temperature $\theta_{e}$ as a function of $M$ and
$\partial_{3}p$ (see Corollary 10). Inserting that expression into the
definition of the potential vorticity $PV_{e}$ produces the first form of
$PV_{e}$-and-$M$ inversion recorded in (1).

Alternatively, hydrostatic balance and the moist variable $M$ may also be used
to write the water $q$ as a function of $M$ and $\partial_{3}p$, indicating
that the sign of $q$ agrees everywhere with the sign of $M-\partial_{3}p.$
This allows to rewrite $PV_{e}$-and-$M$ inversion as in (2), which is the
formulation of $PV_{e}$-and-$M$ inversion used from now on.

Brief summary of PQG literature. References [WSSM19, WSS+20] explain how
$PV_{e}$-and-$M$ inversion can be used to extract the balanced component of a
given atmospheric state, including the balanced parts of water, winds and
temperature. The authors emphasize that the $PV_{e}$-and-$M$ formulation of
PQG is both balanced and invertible, whereas these key ingredients had not
been combined together in previous studies of moist potential vorticity. In
order to achieve balance, the PQG equations were derived formally in [SS17] as
a distinguished asymptotic limit, and the limiting process has since garnered
more attention from the analytical point of view in [ZSS21], and from the
numerical point of view in [ZSS22a, ZSS22b].

Physical phenomena inherent in the PQG system have been explored through
analytical and numerical studies aimed at specific classes of solutions that
are important in meteorological applications [WSS17, ESS20a, WSS19]. When it
comes to the study of the PQG equations “at large”, i.e. away from specific
classes of solutions, only numerical studies have been carried out so far
[ESS20b, HESS21]. It is worth noting that both of the numerical studies
[ESS20b, HESS21] consider the 2-level PQG system, in which the vertical
dependence is simplified to include only two parallel horizontal slices.

##  4 Variational formulation and basic properties of the energy

In this section we begin working towards the proof of Theorems 1.1. We begin
by defining the energy, whose minimizers will be solutions of $PV_{e}$-and-$M$
inversion, precisely. We then study elementary properties of this energy
pertaining to its convexity, coercivity, and differentiability. These
properties will then be used throughout the remainder of this paper. Finally
we conclude this section by proving that indeed minimizers of this energy are
solutions of $PV_{e}$-and-$M$ inversion.

Here is the energy in question.

######  Definition 4.1 (Energy).

For any $M\in L^{2}$, $PV_{e}\in H^{-1}$, and $p\in\mathring{H}^{1}$, where
all spaces are taken over $\mathbb{T}^{3}$ and
$H^{-1}\vcentcolon={(\mathring{H}^{1})}^{*}$, we define

$E(p)\vcentcolon=\int_{\mathbb{T}^{3}}\frac{1}{2}{\lvert\nabla
p\rvert}^{2}-\frac{1}{4}{\min(M-\partial_{3}p,\,0)}^{2}+\langle
PV_{e},\,p\rangle.$ (4)

Here and throughout, unless specified otherwise, the pairing
$\langle\,\cdot\,,\,\cdot\,\rangle$ used in (4) is the
$H^{-1}\times\mathring{H}^{1}$ duality pairing, i.e.
$\langle\Lambda,\,p\rangle=\Lambda(p)$ for every $\Lambda\in H^{-1}$ and
$p\in\mathring{H}^{1}$.

As noted previously, an essential feature of the energy introduced in
Definition 4.1 above is the lack of smoothness due to the presence of the
minimum function. As an intermediate step, it is therefore often useful to
work with a smoothed out version of the problem. This is done by employing
mollifiers. To avoid any confusion with differing conventions we record below
what we mean by a standard mollifier.

######  Definition 4.2 (Mollifier).

A _standard mollifier_ is a non-negative smooth function
$\varphi:\mathbb{R}\to\mathbb{R}$ whose support is contained in $(-1,\,1)$ and
which is _normalised_ , meaning that $\int\varphi=1$. Moreover we say that
$\varphi$ is _centered_ if $\int y\varphi(y)dy=0$ – i.e. the probability
distribution $\varphi(y)dy$ has mean zero. We will often equivalently refer to
$\varphi_{\varepsilon}$ as a mollifier to mean the family of mollifiers
defined by
$\varphi_{\varepsilon}\vcentcolon=\frac{1}{\varepsilon}\varphi\left(\frac{\cdot}{\varepsilon}\right)$.

We now turn our attention to the study of the convexity of the energy. First
we show that the energy density is strongly convex – see Lemma A.1 for a
reminder of equivalent characterisations of strong convexity.

######  Lemma 4.3 (Strong convexity of the energy density).

For any fixed $r\in\mathbb{R}$ the function $e:\mathbb{R}^{3}\to\mathbb{R}$
defined as $e(u)\vcentcolon=\frac{1}{2}{\lvert
u\rvert}^{2}-\frac{1}{4}\min{\left(r-u_{3},\,0\right)}^{2}$ is
$\frac{1}{2}$-strongly convex.

######  Proof 4.4.

Write $e(u)=\frac{1}{2}{\lvert u\rvert}^{2}-\frac{1}{2}f_{0}(r-u_{3})$ for
$f_{0}(s)\vcentcolon=\frac{1}{2}{\min(s,\,0)}^{2}$. For
$\varphi_{\varepsilon}$ a standard mollifier define
$f_{\varepsilon}\vcentcolon=f_{0}*\varphi_{\varepsilon}$ and
$e_{\varepsilon}(u)\vcentcolon=\frac{1}{2}{\lvert
u\rvert}^{2}-\frac{1}{2}f_{\varepsilon}(r-u_{3})$. Then we have that
$\nabla^{2}e_{\varepsilon}(u)=I-\frac{1}{2}f_{\varepsilon}^{\prime\prime}(r-u_{3})e_{3}\otimes
e_{3}$ where ${\lvert\lvert
f_{\varepsilon}^{\prime\prime}\rvert\rvert}_{L^{\infty}}\leqslant\text{Lip}(f_{0}^{\prime})=1$
since $f_{0}^{\prime}(s)=\min(s,\,0)$, and hence
$\nabla^{2}e_{\varepsilon}\geqslant\frac{1}{2}I$ which shows, by virtue of
Lemma A.3 that $e_{\varepsilon}$ is $\frac{1}{2}$-strongly convex. We may then
pass to the limit as $\varepsilon\to 0$ in any of the _first-order_
characterisations of strong convexity of Lemma A.1 to deduce (since $f_{0}\in
C^{1}$, and hence $e\in C^{1}$ such that $e_{\varepsilon}\to g$ in $C^{1}$)
that $e$ itself is $\frac{1}{2}$-strongly convex.

We continue studying the convexity of the energy, now proving that the energy
itself is strongly convex. The weak lower semi-continuity of the energy then
follows. We also record a coercivity estimate for the energy.

######  Proposition 1 (Properties of the energy).

The energy $E$ introduced in Definition 4.1 is $\frac{1}{2}$-strongly convex
and weakly lower semi-continuous over $\mathring{H}^{1}$. Moreover $E$ is
coercive over $\mathring{H}^{1}$ in the sense that
$E(p)\geqslant\frac{1}{16}{\lvert\lvert
p\rvert\rvert}_{\mathring{H}^{1}}^{2}-\frac{3}{4}{\lvert\lvert
M\rvert\rvert}_{L^{2}}^{2}-4{\lvert\lvert PV_{e}\rvert\rvert}_{H^{-1}}^{2}$
for every $p\in\mathring{H}^{1}$.

######  Proof 4.5.

We may write

$E(p)=\int_{\mathbb{T}^{3}}e_{M}\left(x,\,\nabla p(x)\right)dx+\langle
PV_{e},\,p\rangle$ (5)

for $e_{M}(x,\,u)\vcentcolon=\frac{1}{2}{\lvert
u\rvert}^{2}-\frac{1}{4}{\min\left(M(x)-u_{3},\,0\right)}^{2}$. Lemma 4.3 then
tells us that, for any $M\in L^{2}$, $e_{M}(x,\,\cdot\,)$ is
$\frac{1}{2}$-strongly convex for almost every $x$ in $\mathbb{T}^{3}$, and so
the strong convexity of $E$ follows. The weak lower semi-continuity then
follows from Lemma A.6 since strong convexity implies convexity. We now turn
our attention to the coercivity of $E$. We introduce
$H_{u}\vcentcolon=\mathds{1}(M<\partial_{3}p)$ and
$H_{s}\vcentcolon=1-H_{u}=\mathds{1}(M\geqslant\partial_{3}p)$ to rearrange
$E$ such that

$\displaystyle E(p)-\langle
PV_{e},\,p\rangle=\int_{\mathbb{T}^{3}}\frac{1}{2}{\lvert\nabla
p\rvert}^{2}H_{s}+\Bigg{(}\frac{1}{2}{\lvert\nabla_{h}p\rvert}^{2}+\underbrace{\frac{1}{4}{(\partial_{3}p)}^{2}+\frac{1}{2}M\partial_{3}p-\frac{1}{4}M^{2}}_{=\vcentcolon(\star)}\Bigg{)}H_{u}.$
(6)

We note that, by Cauchy-Schwarz,
$(\star)\geqslant\frac{1}{8}{(\partial_{3}p)}^{2}-\frac{3}{4}M^{2}$. Combining
this inequality with (6) above we see that $E(p)-\langle
PV_{e},\,p\rangle\geqslant\frac{1}{8}{\lvert\lvert
p\rvert\rvert}_{\mathring{H}^{1}}^{2}-\frac{3}{4}{\lvert\lvert
M\rvert\rvert}_{L^{2}}$. To conclude we apply Cauchy-Schwarz to $\langle
PV_{e},\,p\rangle$.

######  Remark 4.6 (Non-negativity of the energy).

Proposition 1 shows that the quadratic part of $E$, namely $E-PV_{e}$, may not
necessarily be non-negative. Indeed, for $M=-\partial_{3}p$ and
$\nabla_{h}p=0$ we read from (6) that $E(p)-\langle
PV_{e},\,p\rangle=\int_{\mathbb{T}^{3}}\frac{1}{2}M^{2}(H_{s}-H_{u})$. If one
sought to work with an energy whose quadratic part was non-negative this issue
would be easily remedied by adding $\frac{1}{2}\int_{\mathbb{T}^{3}}M^{2}$ to
the energy: this has no bearing on the minimization of the energy with respect
to $p$ and a Cauchy-Schwarz argument as in the proof of Proposition 1 shows
that then $E-PV_{e}+\frac{1}{2}{\lvert\lvert M\rvert\rvert}_{L^{2}}\geqslant
0$. (This is actually optimal: the example above where $M=-\partial_{3}p$
shows that adding any smaller multiple of the $L^{2}$ norm of $M$ to the
energy would not be sufficient to guarantee the non-negativity of its
quadratic part.)

Having studied the convexity and coercivity of the energy we now turn our
attention towards its differentiability. First we show that the energy is
Gâteaux differentiable and obtain an expression for its derivative.

######  Lemma 4.7 (Gâteaux differentiability of the energy).

The energy $E$ introduced in Definition 4.1 is Gâteaux differentiable on
$\mathring{H}^{1}$ with Gâteaux derivative at any $p\in\mathring{H}^{1}$,
denoted by $DE(p)$, given by

$DE(p)\phi=\int_{\mathbb{T}^{3}}\nabla
p\cdot\nabla\phi+\frac{1}{2}\min(M-\partial_{3}p,\,0)\partial_{3}\phi+\langle
PV_{e},\,\phi\rangle\text{ for every }\phi\in\mathring{H}^{1}.$ (7)

######  Proof 4.8.

Since $PV_{e}\in H^{-1}$ we know that $p\mapsto\langle PV_{e},\,p\rangle$ is
Gâteaux differentiable (and equal to its derivative) and so we focus our
attention on the remaining terms in the energy, namely
$E_{0}(p)\vcentcolon=\int_{\mathbb{T}^{3}}\frac{1}{2}{\lvert\nabla
p\rvert}^{2}-\frac{1}{4}{\min(M-\partial_{3}p,\,0)}^{2}.$ It suffices to show
that $E_{0}$ is Gâteaux differentiable with

$DE_{0}(p)\phi=\int_{\mathbb{T}^{3}}\nabla
p\cdot\nabla\phi+\frac{1}{2}\min(M-\partial_{3}p,\,0)\partial_{3}\phi\text{
for every }\phi\in\mathring{H}^{1}.$ (8)

To do that we fix $p,\,\phi\in\mathring{H}^{1}$ and define

$\tilde{e}(t,\,x)\vcentcolon=\frac{1}{2}{\lvert\nabla
p(x)+t\nabla\phi(x)\rvert}^{2}-\frac{1}{4}{\min\left(M(x)-\left[\partial_{3}p(x)+t\partial_{3}\phi(x)\right],\,0\right)}^{2}$

such that $E_{0}(p+t\phi)=\int_{\mathbb{T}^{3}}\tilde{e}(t,\,x)dx.$ Since
$\tilde{e}(t,\,\cdot\,)\in L^{1}(\mathbb{T}^{3})$ for all $t$ with

$\frac{\partial\tilde{e}}{\partial t}(t,\,x)=\left[\nabla
p(x)+t\nabla\phi(x)\right]\cdot\nabla\phi(x)+\frac{1}{2}\min\left(M(x)-\left[\partial_{3}p(x)+t\partial_{3}\phi(x)\right],\,0\right)\partial_{3}\phi(x)$
(9)

an application of the Dominated Convergence Theorem (more precisely Theorem
2.27 in [Fol99]) tells us that $t\mapsto E_{0}(p+t\phi)$ is differentiable,
with
$\frac{\mathrm{d}}{\mathrm{d}t}E_{0}(t+t\phi)=\int_{\mathbb{T}^{3}}\frac{\partial\tilde{e}}{\partial
t}.$ Evaluating at $t=0$ and using (9) then yields (8).

We now bootstrap up from Gâteaux to Fréchet differentiability in light of the
continuity of the Gâteaux derivative of the energy.

######  Corollary 2 (Fréchet differentiability of the energy).

The Gâteaux derivative of the energy $E$ introduced in Definition 4.1 is
Lipschitz, with specifically the estimate ${\lvert\lvert
DE(p_{1})-DE(p_{2})\rvert\rvert}_{H^{-1}}\leqslant\frac{3}{2}{\lvert\lvert
p_{1}-p_{2}\rvert\rvert}_{\mathring{H}^{1}}$ for every
$p_{1},\,p_{2}\in\mathring{H}^{1}$, and so $E$ is Fréchet differentiable.

######  Proof 4.9.

This follows from the expressions for $DE$ recorded in Lemma 4.7 and a simple
estimate using the fact that the map $s\mapsto\min(s,\,0)$ is $1$-Lipschitz.

With the Fréchet differentiability of the energy in hand we may now obtain a
simple (but very handy!) estimate that shows that as a consequence of its
strong convexity, the energy essentially acts as a norm about its minimizer.

######  Lemma 4.10 (The energy is essentially a norm abouts its minimizer).

Suppose that the energy $E$ introduced in Definition 4.1 has a minimizer
$p^{*}\in\mathring{H}^{1}$. For any $p\in\mathring{H}^{1}$ the following
estimate holds: ${\lvert\lvert
p-p^{*}\rvert\rvert}_{\mathring{H}^{1}}^{2}\leqslant
4\left(E(p)-E(p^{*})\right)=4\left(E(p)-\min E\right).$

######  Proof 4.11.

Corollary 2 tells us that $E$ is Fréchet differentiable and so the result
follows from combining Proposition 1 and Lemma A.1. Since $E$ is
$\frac{1}{2}$-strongly convex and $DE(p^{*})=0$ at the minimizer we obtain
that, for any $p\in\mathring{H}^{1}$,

$E(p)-E(p^{*})\geqslant\langle
DE(p^{*}),\,p-p^{*}\rangle+\frac{1/2}{2}{\lvert\lvert
p-p^{*}\rvert\rvert}_{\mathring{H}^{1}}^{2}=\frac{1}{4}{\lvert\lvert
p-p^{*}\rvert\rvert}_{\mathring{H}^{1}}^{2},$

as desired.

We now conclude this section by proving what is perhaps its main result,
namely that minimizers of the energy introduced in Definition 4.1 are
precisely the solutions of $PV_{e}$-and-$M$ inversion.

######  Lemma 4.12 (Variational formulation of $PV_{e}$-and-$M$ inversion).

For any $p\in\mathring{H}^{1}$, $p$ is a global minimizer of $E$ introduced in
Definition 4.1 if and only if it is an $\mathring{H}^{1}$-weak solution of

$-\Delta p-\frac{1}{2}\partial_{3}\left(\min(M-\partial_{3}p,\,\
0)\right)=-PV_{e}$ (10)

in the sense that

$\int_{\mathbb{T}^{3}}\nabla
p\cdot\nabla\phi+\frac{1}{2}\min(M-\partial_{3}p,\,0)(\partial_{3}\phi)=-\langle
PV_{e},\,\phi\rangle\text{ for every }\phi\in\mathring{H}^{1}.$ (11)

######  Proof 4.13.

Lemma 4.7 tells us that (11) is equivalent to $DE(p)=0$. Since $E$ is
differentiable and strictly convex (c.f. Corollary 2 and Proposition 1) the
claim then follows immediately.

##  5 Well-posedness

In this section we use the variational formulation of $PV_{e}$-and-$M$
inversion to show that the problem is well-posed. We first show that unique
solutions exist in Theorem 5.1 and then turn our attention to their
regularity.

As discussed previously, the de Giorgi theory is central to that matter of
regularity and so the first main result of this section is Theorem 5.7 in
which the Hölder regularity of the gradient of solutions of $PV_{e}$-and-$M$
inversion is established. With this key result in hand we then turn our
attention to the matter of higher regularity, establishing sufficient
conditions on $PV_{e}$ and $M$ for the solutions to be classical, in each
phase, in Theorem 5.12 and finally establishing sufficient conditions on
$PV_{e}$ and $M$ for the solutions to be smooth in each phase in Theorem 5.14.

We conclude this section with a discussion of how to construct explicit one-
dimensional ($x_{3}$–dependent) solutions, which in particular show that the
regularity theory developed here is sharp.

First, as promised, we establish existence and uniqueness.

######  Theorem 5.1 (Existence and uniqueness).

For every $M\in L^{2}$, and $PV_{e}\in H^{-1}$ the energy $E$ introduced in
Definition 4.1 has a unique global minimizer in $\mathring{H}^{1}$ which is
also the unique weak solution of (11).

######  Proof 5.2.

For any $p\in\mathring{H}^{1}$, $M\in L^{2}$, and $PV_{e}\in H^{-1}$ the
energy $E(p)$ is finite and the coercivity of the energy recorded in
Proposition 1 tells us that $E(p)$ is bounded below independently of $p$.
Therefore $\inf_{\mathring{H}^{1}}E$ is well-defined. We may thus take a
minimizing sequence $(p_{n})\subseteq\mathring{H}^{1}$ satisfying
$E(p_{n})\downarrow\inf E$ as $n\to\infty$. The coercivity of $E$ implies that
the sequence $(p_{n})$ is bounded in $\mathring{H}^{1}$, and so we deduce from
Kakutani’s theorem that there is a subsequence $(p_{n_{k}})$ which is
$\mathring{H}^{1}$-weakly convergent to some $p_{\infty}\in\mathring{H^{1}}$.
Finally, since Proposition 1 tells us that the energy is weakly lower semi-
continuous in $\mathring{H}^{1}$, we deduce that
$E(p_{\infty})\leqslant\liminf_{k\to\infty}E(p_{n_{k}})=\inf E,$ i.e.
$p_{\infty}$ is a minimizer of $E$, as desired. In particular the strict
convexity of $E$ established in Proposition 1 guarantees that this minimizer
is unique. Lemma 4.12 then ensures that this unique minimizer is precisely the
unique weak solution of (11).

We now turn our attention to the regularity of solutions of $PV_{e}$-and-$M$
inversion. First we obtain that if $PV_{e}$ and $M$ are sufficiently regular
then the solutions lie in $H^{2}$. This argument relies on the use of finite
differences (see Definition B.1 and Lemma B.2).

######  Lemma 5.3 ($H^{2}$ regularity).

Let $M\in H^{1}$ and $PV_{e}\in L^{2}$. Any $\mathring{H}^{1}$-weak solution
$p$ of (11) belongs to $H^{2}$ and satisfies ${\lvert\lvert
p\rvert\rvert}_{\dot{H}^{2}}\leqslant{\lvert\lvert
M\rvert\rvert}_{\dot{H}^{1}}+2{\lvert\lvert PV_{e}\rvert\rvert}_{L^{2}}.$

######  Proof 5.4.

We define, for $g_{0}(s)\vcentcolon=\min(s,\,0)$, and for any
$x\in\mathbb{T}^{3}$, the function
$\zeta_{h}^{i}(x)\vcentcolon=\int_{0}^{1}g_{0}^{\prime}\left((I+\theta\Delta_{h}^{i})(M-\partial_{3}p)(x)\right)d\theta$
where recall that the finite difference operator $\Delta_{h}^{i}$ is
introduced in Definition B.1, and note that this function is well-defined
almost everywhere independently of the representative of $M-\partial_{3}p\in
L^{2}$ used since for any fixed $h$ and $i$, defining $\zeta_{h}^{i}$ only
requires the evaluation of $M-\partial_{3}p$ and
$(M-\partial_{3}p)(\,\cdot\,+he_{i})$. Note that $\zeta$ is defined to be
precisely the term that will appear, by virtue of the chain rule for finite
differences recorded in Lemma B.2, when applying a finite difference to the
function $\min(M~{}-~{}\partial_{3}p,\,0)~{}=~{}g_{0}(M~{}-~{}\partial_{3}p).$
Note also that, since the Lipschitz constant of $g_{0}$ is equal to one,
$\lvert\zeta_{h}^{i}(x)\rvert\leqslant 1$ for every $x\in\mathbb{T}^{3}$.

Now we use $\Delta_{-h}^{i}\Delta_{h}^{i}p$ as a test function (i.e.
${(\Delta_{h}^{i})}^{*}\Delta_{h}^{i}$, as indicated by the computation of the
adjoint of the finite difference operator in Lemma B.2). This produces, using
Lemma B.2,

$\int(PV_{e})\Delta_{-h}^{i}\Delta_{h}^{i}p=\int{\lvert\nabla\Delta_{h}^{i}p\rvert}^{2}-\frac{1}{2}{(\partial_{3}\Delta_{h}^{i}p)}^{2}\zeta_{h}^{i}+\frac{1}{2}(\Delta_{h}^{i}M)(\partial_{3}\Delta_{h}^{i}p)\zeta_{h}^{i}.$

Therefore, since $\lvert\zeta_{h}^{i}\rvert\leqslant 1$,

$\displaystyle\frac{1}{2}\int{\lvert\nabla\Delta_{h}^{i}p\rvert}^{2}$
$\displaystyle\leqslant\int{\lvert\nabla\Delta_{h}^{i}p\rvert}^{2}-\frac{1}{2}{(\partial_{3}\Delta_{h}^{i}p)}^{2}\zeta_{h}^{i}$
$\displaystyle=\int-\frac{1}{2}(\Delta_{h}^{i}M)(\partial_{3}\Delta_{h}^{i}p)\zeta_{h}^{i}+\int(\Delta_{h}^{i}PV_{e})(\Delta_{h}^{i}p)$
$\displaystyle\leqslant\left(\frac{1}{2}{\lvert\lvert\Delta_{h}^{i}M\rvert\rvert}_{L^{2}}+{\lvert\lvert\Delta_{h}^{i}PV_{e}\rvert\rvert}_{H^{-1}}\right){\lvert\lvert\nabla\Delta_{h}^{i}p\rvert\rvert}_{L^{2}}$

such that, using Lemma B.2 once more,

${\lvert\lvert\Delta_{h}^{i}\nabla
p\rvert\rvert}_{L^{2}}\leqslant{\lvert\lvert\Delta_{h}^{i}M\rvert\rvert}_{L^{2}}+2{\lvert\lvert\Delta_{h}^{i}PV_{e}\rvert\rvert}_{H^{-1}}\leqslant\lvert
h\rvert\left({\lvert\lvert\nabla M\rvert\rvert}_{L^{2}}+2{\lvert\lvert
PV_{e}\rvert\rvert}_{L^{2}}\right).$

So finally we use Lemma B.2 one last time to conclude that $\nabla^{2}p\in
L^{2}$, with the estimate ${\lvert\lvert
p\rvert\rvert}_{\dot{H}^{2}}={\lvert\lvert\nabla^{2}p\rvert\rvert}_{L^{2}}\leqslant{\lvert\lvert
M\rvert\rvert}_{\dot{H}^{1}}+2{\lvert\lvert PV_{e}\rvert\rvert}_{L^{2}}.$

The next step in climbing the regularity ladder is to prove that the gradient
of solutions to $PV_{e}$-and-$M$ inversion is Hölder continuous. Following de
Giorgi, we do this by considering the equation satisfied by the gradient of
such solutions. That equation is recorded in the result below.

######  Lemma 5.5 (Equation satisfied by the gradient).

Let $M\in H^{1}$ and $PV_{e}\in L^{2}$. If $p\in\mathring{H}^{1}$ is an
$\mathring{H}^{1}$-weak solution of (10) then $u_{i}\vcentcolon=\partial_{i}p$
is an $\mathring{H}^{1}$-weak solution of

$-\Delta_{h}u_{i}-\frac{1}{2}\partial_{3}\Big{(}[1+\mathds{1}(M\geqslant\partial_{3}p)](\partial_{3}u_{i})\Big{)}=-\partial_{i}PV_{e}+\frac{1}{2}\partial_{3}\Big{(}\mathds{1}(M<\partial_{3}p)(\partial_{i}M)\Big{)}$
(12)

in the sense that, for every $\phi\in\mathring{H}^{1}$,

$\displaystyle\int_{\mathbb{T}^{3}}\nabla_{h}u_{i}\cdot\nabla_{h}\phi+\frac{1}{2}[1+\mathds{1}(M\geqslant\partial_{3}p)](\partial_{3}u_{i})(\partial_{3}\phi)$
$\displaystyle=-\langle\partial_{i}PV_{e},\,\phi\rangle-\int_{\mathbb{T}^{3}}\frac{1}{2}\mathds{1}(M<\partial_{3}p)(\partial_{i}M)(\partial_{3}\phi).$

######  Proof 5.6.

As usual we write $g_{0}(s)\vcentcolon=\min(0,\,s)$. Since $M\in H^{1}$ and
$PV_{e}\in L^{2}$ Lemma 5.3 tells us that $p\in H^{2}$, and so
$M-\partial_{3}p\in H^{1}$. Since $g_{0}$ is Lipschitz and $g_{0}(0)=0$ we
deduce that $g_{0}(M-\partial_{3}p)$ belongs to $H^{1}$ (see Theorem 12.69 of
[Leo17]), with
$\partial_{i}\left(g_{0}(M-\partial_{3}p)\right)=g_{0}^{\prime}(M-\partial_{3}p)(\partial_{i}M-\partial_{3}\partial_{i}p)=\mathds{1}(M<\partial_{3}p)(\partial_{i}M-\partial_{3}\partial_{i}p).$
In particular, for any $\psi\in\mathring{C}^{\infty}$ (the space of smooth
functions with average zero) we may use $\partial_{i}\psi$ as a test function
to see that, using the weak formulation recorded in Lemma 4.12 and integrating
by parts,

$\displaystyle-\langle\partial_{i}PV_{e},\,\psi\rangle={\left(PV_{e},\,\partial_{i}\psi\right)}_{L^{2}}=\int_{\mathbb{T}^{3}}\nabla\partial_{i}p\cdot\nabla\psi+\frac{1}{2}g_{0}^{\prime}(M-\partial_{3}p)(\partial_{i}M-\partial_{3}\partial_{i}p)(\partial_{3}\psi)$
$\displaystyle=\int_{\mathbb{T}^{3}}\nabla_{h}u_{i}\cdot\nabla_{h}\psi+\left(1-\frac{1}{2}\mathds{1}(M<\partial_{3}p)\right)(\partial_{3}u_{i})(\partial_{3}\psi)+\frac{1}{2}\mathds{1}(M<\partial_{3}p)(\partial_{i}M)(\partial_{3}\psi)$

and so indeed

$\displaystyle\int_{\mathbb{T}^{3}}\nabla_{h}u_{i}\cdot\nabla_{h}\psi+\frac{1}{2}\left(1+\mathds{1}(M\geqslant\partial_{3}p)\right)(\partial_{3}u_{i})(\partial_{3}\psi)$
$\displaystyle=\int_{\mathbb{T}^{3}}-\frac{1}{2}\mathds{1}(M<\partial_{3}p)(\partial_{i}M)(\partial_{3}\psi)-\langle\partial_{i}PV_{e},\,\psi\rangle.$
(13)

In particular, since $u_{i}=\partial_{i}p\in H^{1}$, $\partial_{i}M\in L^{2}$,
and $\partial_{i}PV_{e}\in H^{-1}$ we may for any $\phi\in\mathring{H}^{1}$
find an approximating sequence
${(\psi_{n})}_{n}\subseteq\mathring{C}^{\infty}$ converging to $\psi$ in
$\mathring{H}^{1}$ such that (13) holds for $\psi=\psi_{n}$ and hence, passing
to the limit, (13) also holds for $\psi=\phi$, as desired.

We are now equipped to prove the first main result of this section, which
states that the gradient of solutions of $PV_{e}$-and-$M$ inversion is Hölder
continuous provided that $PV_{e}$ and $M$ are sufficiently regular. The proof
of this result makes uses of classical ideas in the sense that it relies on
framing the de Giorgi argument in a manner compatible with the Campanato
formulation of Hölder continuity (see Definition B.4 and Lemma B.6 ). We
include the proof here for the following reason: a plethora of expository work
(of which we only mention [CV10, Vas16]) discusses the _homogeneous_ de Giorgi
method in the elliptic case where there is no forcing term and some textbooks
such [HL11] discuss the _inhomogeneous_ case by analogy with the classical
treatment of elliptic equations with forcing, but an explicit proof, both
sufficiently detailed and simple for the non-specialist, is difficult to find.
The intent is that the proof provided here fits that bill. (Note that
auxiliary results used in the proof are relegated to Appendix B – we emphasize
the assembly of these ideas here.)

######  Theorem 5.7 (Hölder continuity of the gradient).

Suppose that $M\in H^{1}$ satisfies $\nabla M\in
L^{q}(\left\\{M<\partial_{3}p\right\\})$ and that $PV_{e}\in L^{q}$ for some
$q>d=3$. There exists $\alpha=\alpha(q)\in(0,\,1)$ with $\alpha\leqslant
1-\frac{d}{q}$ such that if $p$ is an $\mathring{H}^{1}$-weak solution of (10)
then $\nabla p\in C^{0,\,\alpha}$ and the following estimate holds:

${\lvert\lvert\nabla p\rvert\rvert}_{C^{0,\,\alpha}(\mathbb{T}^{3})}\leqslant
C\left({\lvert\lvert\nabla
M\rvert\rvert}_{L^{q}(M<\partial_{3}p)}+{\lvert\lvert
PV_{e}\rvert\rvert}_{L^{q}(\mathbb{T}^{3})}\right)$ (14)

where $C=C(q)>0$.

######  Proof 5.8.

Fix $i=1,\,2,\,3$ and write $u\vcentcolon=\partial_{i}p$. Note that, as
derived in Lemma 5.5, $u$ is an $\mathring{H}^{1}$-weak solution of (12) which
may be written as

$-\nabla\cdot(A(x)\nabla u)=\nabla\cdot F$ (15)

for
$A\vcentcolon=\widetilde{I}+\frac{1}{2}\left(1+\mathds{1}(\partial_{3}p<M)\right)e_{3}\otimes
e_{3}$ and
$F\vcentcolon=PV_{e}e_{i}+\frac{1}{2}\mathds{1}(M<\partial_{3}p)(\partial_{i}M)e_{3},$
where $\widetilde{I}=I-e_{3}\otimes e_{3}$. Crucially: $A\in L^{\infty}$ is
uniformly elliptic since
$A(x)\xi\cdot\xi\geqslant\frac{1}{2}{\lvert\xi\rvert}^{2}$ for all
$x\in\mathbb{R}^{3}$.

We will now use (15) to prove that, for any $x_{0}\in\mathbb{T}^{3}$ and
$r>0$,

$\int_{B(x_{0},\,r)}{\lvert\nabla u\rvert}^{2}\leqslant C{\lvert\lvert
F\rvert\rvert}_{L^{q}}^{2}r^{d-2+2\alpha}$ (16)

for some $\alpha=\alpha(q)\in(0,\,1)$ and $C=C(q)>0$ (recall that $d=3$). It
will then follow from Proposition 6 that $u\in C^{0,\,\alpha}$, with the
desired estimate (14).

So let us fix $x_{0}\in\mathbb{T}^{3}$ and $r>0$. To establish (16) we split
the solution $u$, inside the ball $B=B(x_{0},\,r)$, into two parts. One part,
$v$, accounts for the forcing term $\nabla\cdot F$ while the other part, $w$,
solves a homogeneous problem (i.e. without forcing) in order to correct for
the homogeneous Dirichlet boundary conditions imposed on $v$. More precisely
we define $v$ to be the $\mathring{H}^{1}$-weak solution of

$\left\\{\begin{aligned} &-\nabla\cdot(A(x)\nabla v)=\nabla\cdot F&\text{ in
}B\text{ and }\\\ &v=0&\text{ on }\partial B,\end{aligned}\right.$

i.e. $v\in H_{0}^{1}(B)$, and define $w\vcentcolon=u-v$, such that indeed
$u=v+w$, which means that $w$ is an $H^{1}$-weak solution of
$-\nabla\cdot(A(x)\nabla w)=0$ in $B$ (but we know nothing, and need to know
nothing, about the boundary values of $w$). Combining Propositions 8 and 9 we
deduce that there exist $\alpha_{1}$ and $\alpha_{2}=1-\frac{d}{q}$ both in
$(0,\,1)$, such that, for any $0<\rho\leqslant r$,

$\displaystyle\int_{B(x_{0},\,\rho)}{\lvert\nabla u\rvert}^{2}$
$\displaystyle\leqslant 2\left(\int_{B(x_{0},\,\rho)}{\lvert\nabla
w\rvert}^{2}+\int_{B(x_{0},\,\rho)}{\lvert\nabla v\rvert}^{2}\right)$
$\displaystyle\leqslant
C\left({\left(\frac{\rho}{r}\right)}^{d-2+2\alpha_{1}}\int_{B(x_{0},\,r)}{\lvert\nabla
w\rvert}^{2}+\int_{B(x_{0},\,r)}{\lvert\nabla v\rvert}^{2}\right)$
$\displaystyle\leqslant
C\left({\left(\frac{\rho}{r}\right)}^{d-2+2\alpha_{1}}\int_{B(x_{0},\,r)}{\lvert\nabla
u\rvert}^{2}+r^{d-2+2\alpha_{2}}{\lvert\lvert
F\rvert\rvert}_{L^{q}}^{2}\right).$

We are now in a position to make use of Lemma B.21 with
$\phi(\rho)=\int_{B(x_{0},\,\rho)}{\lvert\nabla u\rvert}^{2}$,
$B=C{\lvert\lvert F\rvert\rvert}_{L^{q}}^{2}$,
$\beta=d-2+2\min(\alpha_{1},\,\alpha_{2})$, and $\gamma=d-2+2\alpha_{1}$. We
deduce that there exists $\delta\in(\beta,\,\gamma)$ such that, for
$\alpha=\min(\alpha_{1},\,\alpha_{2})$, and since $0<\rho\leqslant R_{0}<1$,

$\displaystyle\int_{B(x_{0},\,\rho)}{\lvert\nabla u\rvert}^{2}\leqslant
C\left[{\left(\frac{\rho}{R_{0}/2}\right)}^{\delta}\phi\left(\frac{R_{0}}{2}\right)+B\rho^{\beta}\right]\leqslant
C\left[\phi\left(\frac{R_{0}}{2}\right)+B\right]\rho^{\beta}$
$\displaystyle\leqslant C\left(\int_{\mathbb{T}^{3}}{\lvert\nabla
u\rvert}^{2}+{\lvert\lvert
F\rvert\rvert}_{L^{q}}^{2}\right)\rho^{d-2+2\alpha}.$ (17)

In particular note that $\alpha\leqslant\alpha_{2}=1-\frac{d}{q}$. To conclude
it suffices to use Lemma B.23 and Hölder’s inequality, such that
$\int_{\mathbb{T}^{3}}{\lvert\nabla u\rvert}^{2}\leqslant C(q){\lvert\lvert
F\rvert\rvert}_{L^{q}}^{2},$ and hence obtain (16) from (17). Proposition 6
then tells us that ${\lvert\lvert u\rvert\rvert}_{C^{0,\,\alpha}}\leqslant
C(q)({\lvert\lvert F\rvert\rvert}_{L^{q}}+{\lvert\lvert
u\rvert\rvert}_{L^{2}})$ and so we may use Lemma B.23 once again to finally
obtain (14).

As is often the case in elliptic regularity theory, having established that
the gradient of solutions is Hölder continuous we immediately deduce useful
consequences and higher regularity. First, a useful consequence: if $PV_{e}$
and $M$ are sufficiently regular then both the unsaturated phase
$\left\\{M<\partial_{3}p\right\\}$ and the saturated phase
$\left\\{\partial_{3}p<M\right\\}$ are open sets.

######  Corollary 3 (Each phase is an open set).

Let $M\in H^{1}$ be continuous with $\nabla M\in L^{q}(M<\partial_{3}p)$ and
let $PV_{e}\in L^{q}$ for some $q>d=3$. If $p$ is an $\mathring{H}^{1}$-weak
solution of (10) then $\partial_{3}p-M$ is continuous and so both the
unsaturated phase $\left\\{M<\partial_{3}p\right\\}$ and the saturated phase
$\left\\{\partial_{3}p<M\right\\}$ are open sets.

######  Proof 5.9.

This follows immediately from Theorem 5.7 which tells us that, under these
assumptions on $M$ and $PV_{e}$, $\nabla p$ is Hölder continuous.

We now continue deducing further results from the Hölder continuity of the
gradient obtained in Theorem 5.7 above. Here we obtain higher-order
regularity, namely sufficient conditions for the Hölder continuity of the
_Hessian_ of solutions of $PV_{e}$-and-$M$ inversion.

######  Lemma 5.10 (Hölder continuity of the Hessian away from the interface).

Let $M\in H^{1}$ be continuous with both $\nabla M$ and $\partial_{3}\nabla M$
belonging to $L^{q}\left(M<\partial_{3}p\right)$, let $PV_{e}\in L^{q}$, and
let $\nabla PV_{e}\in L^{q}\left(M\neq\partial_{3}p\right)$ for some $q>d=3$.
If $p$ is an $\mathring{H}^{1}$-weak solution of (10) then $\nabla^{2}p\in
C^{0,\,\alpha}_{\text{loc}}\left(M\neq\partial_{3}p\right)$ for
$\alpha=1-\frac{d}{q}$.

######  Proof 5.11.

Lemma 5.5 tells us that $u=\partial_{i}p$ is an $\mathring{H}^{1}$-weak
solution of (12). Since Corollary 3 tells us that
$\left\\{M<\partial_{3}p\right\\}$ is an open set we may test (12) against any
$\varphi\in C_{c}^{\infty}(M<\partial_{3}p)$ to obtain

$\int_{\left\\{M<\partial_{3}p\right\\}}\nabla_{h}u\cdot\nabla_{h}\varphi+\frac{1}{2}(\partial_{3}u)(\partial_{3}\varphi)=-\langle\partial_{i}PV_{e},\,\varphi\rangle-\int_{\left\\{M<\partial_{3}p\right\\}}\frac{1}{2}(\partial_{i}M)(\partial_{3}\varphi)$

such that, by passing to the limit from $C_{c}^{\infty}$ to $H_{0}^{1}$, we
obtain that $u\in H^{1}$ is an $H_{0}^{1}$-weak solution of

$\left(-\Delta_{h}-\frac{1}{2}\partial_{3}^{2}\right)u=-\partial_{i}PV_{e}+\frac{1}{2}\partial_{3}\partial_{i}M\text{
in }\left\\{M<\partial_{3}p\right\\}$

(note that $u$ does _not_ belong to $H_{0}^{1}$, it is only the test functions
that belong to $H_{0}^{1}$). By classical theory (the elliptic operator has
_constant coefficients_ when restricted to a single phase, so we may use for
example Theorem 3.13 of [HL11]) it follows that, since
$\frac{1}{2}\partial_{3}\partial_{i}M-\partial_{i}PV_{e}\in
L^{q}\left(M<\partial_{3}p\right)$, consequently $\nabla u\in
C^{0,\,\alpha}_{\text{loc}}\left(M<\partial_{3}p\right)$ for
$\alpha=1-\frac{d}{q}$. Similarly we obtain that $u\in H^{1}$ is an
$H_{0}^{1}$-weak solution of

$-\Delta u=-\partial_{i}PV_{e}\text{ in }\left\\{\partial_{3}p<M\right\\}$

such that, since $\partial_{i}PV_{e}\in L^{q}\left(\partial_{3}p<M\right)$, it
follows that $\nabla u\in
C^{0,\,\alpha}_{\text{loc}}\left(\partial_{3}p<M\right)$ as desired.

We are now ready to prove the second main result of this section, namely
identifying sufficient conditions on $PV_{e}$ and $M$ which ensure that
solutions of $PV_{e}$-and-$M$ inversion are classical in each phase.

######  Theorem 5.12 (Classical solutions in each phase).

Let $M\in H^{1}$ be continuous with both $\nabla M$ and $\nabla\partial_{3}M$
belonging to $L^{q}\left(M<\partial_{3}p\right)$, let $PV_{e}\in L^{q}$, and
also let $\nabla PV_{e}\in L^{q}\left(M\neq\partial_{3}p\right)$ for some
$q>d=3$. If $p$ is an $\mathring{H}^{1}$-weak solution of (10) then it is a
classical solution of

$\displaystyle\left(-\Delta_{h}-\frac{1}{2}\partial_{3}^{2}\right)p$
$\displaystyle=-PV_{e}+\frac{1}{2}\partial_{3}M\text{ in
}\left\\{M<\partial_{3}p\right\\}$ (18) and $\displaystyle-\Delta p$
$\displaystyle=-PV_{e}\text{ in }\left\\{\partial_{3}p<M\right\\}.$ (19)

######  Proof 5.13.

Corollary 3 tells us that $\left\\{M<\partial_{3}p\right\\}$ is open so we may
test (10) against $\varphi\in C_{c}^{\infty}(M<\partial_{3}p)$, obtaining

$\displaystyle-\int_{\left\\{M<\partial_{3}p\right\\}}(PV_{e})\varphi=\int_{\left\\{M<\partial_{3}p\right\\}}\nabla_{h}p\cdot\nabla_{h}\varphi+\frac{1}{2}(\partial_{3}p)(\partial_{3}\varphi)+\frac{1}{2}M(\partial_{3}\varphi).$
(20)

The usual approximation argument then shows that (20) actually holds for all
$\varphi\in H_{0}^{1}(M<\partial_{3}p)$. Therefore, since Lemma 5.10 tells us
that $p\in C^{2}(M<\partial_{3}p)$ while Morrey’s embedding tells us that, for
$q>d$, $PV_{e},\,\partial_{3}M\in C^{0,\,1-d/q}_{\text{loc}}(M<\partial_{3}p)$
we conclude that indeed $p$ is a classical solution of (18). In the same way
we deduce from (10) and $PV_{e}$ having average zero that

$\int_{\left\\{\partial_{3}p<M\right\\}}\nabla
p\cdot\nabla\phi+(PV_{e})\phi=0\text{ for all }\phi\in
H_{0}^{1}(\partial_{3}p<M),$

from which it follows as above that $p$ is a classical solution of (19).

We now record the penultimate main result of this section, which consists in a
higher-regularity result as well as in establishing that, if both $PV_{e}$ and
$M$ are smooth, then so is, in each phase, the solution of $PV_{e}$-and-$M$
inversion.

######  Theorem 5.14 (Higher regularity in each phase).

Let $q>d=3$, let $k\geqslant 0$ be an integer, and let $\alpha\in(0,\,1)$.
Suppose that

  * •

$M\in H^{1}$ with $\nabla M,\,\nabla\partial_{3}M\in L^{q}(M<\partial_{3}p)$
and $\partial_{3}M\in C^{k,\,\alpha}(M<\partial_{3}p)$ and that

  * •

$PV_{e}\in L^{q}(\mathbb{T}^{3})$ with $\nabla PV_{e}\in
L^{q}(M\neq\partial_{3}p)$ and $PV_{e}\in C^{k,\,\alpha}(M\neq\partial_{3}p)$.

If $p$ is an $\mathring{H}^{1}$-weak solution of (10) then $p\in
C^{k,\,\alpha}_{\text{loc}}(M\neq\partial_{3}p)$. In particular if
$\partial_{3}M$ and $PV_{e}$ are $C^{\infty}$ in
$\left\\{M<\partial_{3}p\right\\}$ and $\left\\{M\neq\partial_{3}p\right\\}$,
respectively, then $p$ is $C^{\infty}$ in
$\left\\{M\neq\partial_{3}p\right\\}$.

######  Proof 5.15.

Under these assumptions Theorem 5.12 tells us that $p$ is a classical solution
of (18)–(19). The result then follows from classical Schauder interior
estimates (see for example Theorem 6.17 of [GT01]).

As we head towards the conclusion of this section we provide a recipe for the
construction of explicit one-dimensional ($x_{3}$–dependent) solutions. This
construction is of particular interest since it can then be used to find
examples where $PV_{e}$ and $M$ are both smooth and yet the solution of
$PV_{e}$-and-$M$ inversion is not. More precisely, that solution will be
$C^{1,\,1}$ but not $C^{2}$, which shows that Theorem 5.7 above is sharp.

Two remarks are in order. First, note that below we specify that the period of
$PV_{e}$ and $M$ is $2\pi$. This is done solely to fix notation and the same
construction carries through for any period. Second, note that we are actually
only able to specify the _profile_ of $M$. Indeed, the solution constructed
satisfies $PV_{e}$-and-$M$ inversion where the data is $PV_{e}$ and $M-c$,
shifted by some constant $c$ which depends on $PV_{e}$ and $M$.

######  Lemma 5.16 (Explicit one-dimensional solutions).

Suppose that $PV_{e}=PV_{e}(x_{3})$ and $M=M(x_{3})$ are smooth
$2\pi$–periodic functions. Let
$\phi_{m}(x)\vcentcolon=x+\frac{1}{2}\min_{0}(m-x)$ where, as usual,
$\min_{0}(x)\vcentcolon=\min(x,\,0)$. $\phi_{m}$ is invertible for any
$m\in\mathbb{R}$ and the following holds. Let $A$ be an antiderivative of
$PV_{e}$ such that $A^{\prime}=PV_{e}$, let
$\widetilde{\Theta}\vcentcolon=\phi_{M}^{-1}\circ A$, and let
$c\vcentcolon=\fint_{-\pi}^{\pi}\widetilde{\Theta}$. Then
$p=p(x_{3})\vcentcolon=\int_{-\pi}^{x_{3}}(\widetilde{\Theta}-c)$ solves (10)
with data $PV_{e}$ and $M-c$.

######  Proof 5.17.

That $\phi_{m}$ is invertible follows directly from the fact that it is
strictly increasing (since $\phi_{m}^{\prime}\geqslant 1/2>0$). Nonetheless it
will be useful to have an explicit representation of its inverse. A direct
computation shows that its inverse is given by
$\phi_{m}^{-1}(x)=x-\min_{0}(m-x)$. We may then compute, since
$\phi_{m-c}(x-c)=\phi_{m}(x)-c$, that the identity
$p^{\prime}+\frac{1}{2}{\min}_{0}\,\left((M-c)-p^{\prime}\right)=\phi_{M-c}(p^{\prime})=(\phi_{M}\circ\widetilde{\Theta})-c$
holds. Therefore

$\displaystyle\Delta
p+\frac{1}{2}\partial_{3}\left({\min}_{0}\,((M-c)-\partial_{3}p)\right)=\left[p^{\prime}+\frac{1}{2}{\min}_{0}\,\left((M-c)-p^{\prime}\right)\right]^{\prime}$
$\displaystyle=(\phi_{M}\circ\widetilde{\Theta})^{\prime}=A^{\prime}=PV_{e},$

as desired.

As an immediate consequence we may conclude this section with its last main
result, verifying that the global regularity (including across the interface)
deduced in Theorem 5.7 is sharp. Once again we work with $2\pi$–periodic
functions to fix notation but the same result holds for arbitrary periods.

######  Corollary 4 (Sharpness of the regularity theory).

Consider a smooth periodic function $M=M(x_{3}):(-\pi,\,\pi)\to\mathbb{R}$
such that $M^{\prime}(0)\neq 0$ and
$\operatorname*{sign}M(x)=\operatorname*{sign}x$ for every $x\in(-\pi,\,\pi)$.
In particular this means that $M(0)=0$ (and that $M^{\prime}(0)>0$). We
introduce $c\vcentcolon=\fint_{-\pi}^{\pi}\min_{0}M$ and define
$\Theta\vcentcolon=c-\min_{0}M$. Then the function
$p=p(x_{3})\vcentcolon=\int_{-\pi}^{x_{3}}\Theta$ is a solution of (10) with
data $PV_{e}\equiv 0$ and $M+c$. Moreover the gradient of $p$ is Lipschitz
continuous everywhere but fails to be differentiable at zero.

######  Proof 5.18.

That $p$ is a solution of $PV_{e}$-and-$M$ inversion follows immediately from
Lemma 5.16 since $-\min_{0}M=\phi_{M}^{-1}(0)$ (for $\phi$ as in Lemma 5.16).
Moreover, since $\operatorname*{sign}M=\operatorname*{sign}x$, we have that

$p^{\prime}(x)=\Theta(x)=c-\left\\{\begin{aligned} &M(x)&&\text{if }x<0\text{
and }\\\ &0&&\text{if }x\geqslant 0.\end{aligned}\right.$

We may then conclude that indeed $p^{\prime}$ is Lipschitz continuous since
$M(0)=0$, and yet fails to be differentiable at zero since $M^{\prime}(0)\neq
0$.

A simple example of Corollary 4 in action is provided below.

######  Example 5.19.

Consider $M:(-\pi,\,\pi)\to\mathbb{R}$ given by $M(x_{3})=\sin(x_{3})-1/\pi$.
Then, for $PV_{e}\equiv 0$, the function $p=p(x_{3})$ given by

$p(x_{3})=-\frac{x}{\pi}+\left\\{\begin{aligned} &\cos x&&\text{if }x\leqslant
0\text{ and }\\\ &1&&\text{if }x>0\end{aligned}\right.$

is a solution of $PV_{e}$-and-$M$ inversion (10). Indeed:
$p^{\prime}(x_{3})=-\min(\sin x_{3},\,0)-1/\pi$ and hence
$\left[p^{\prime}+\frac{1}{2}\min(M-p^{\prime},\,0)\right]^{\prime}=0,$ as
desired. In particular we see that $p$ is $C^{1,\,1}$ but not $C^{2}$. This
so-called “baseball cap” example is depicted in Figure 2.

Note that this example is constructed by following the recipe provided above
in Lemma 5.16, using $\sin x_{3}$ as profile for $M$ such that the constant
$c$ is given by $c=-\fint_{-\pi}^{\pi}\min_{0}(\sin x_{3})\,dx_{3}=1/\pi$.

Figure 2:  The “baseball cap”: a very simple example where the gradient of $p$
(i.e. $\Theta$) is Lipschitz continuous but not differentiable:
$M=\sin(x_{3})-1/\pi$ and $PV_{e}=0$. The name refers to the profile of the
equivalent potential temperature $\Theta$.

Note that with Corollary 4 in hand we have now proved every part of Theorem
1.1, the main result of this paper.

######  Proof 5.20 (Proof of Theorem 1.1).

We combine Theorems 5.1, 5.7, and 5.14 and Corollaries 3 and 4.

##  6 Smoothing

The purpose of this section is two-fold. First it details how using
$PV_{e}$-and-$M$ inversion in its form (1) which explicitly involves
Heavisides makes it difficult to determine whether that inversion is an
elliptic problem or not. This is done by showing that so-called agnostic
smoothing of the Heavisides, which does not take into account their dependence
on the moist variable $M$ and the vertical gradient of the pressure
$\partial_{3}p$, fails to determine whether or not $PV_{e}$-and-$M$ inversion
is elliptic.

Second we introduce in this section a more careful way to smooth out the
Heavisides which relies on the variational formulation of $PV_{e}$-and-$M$
inversion. The resulting regularised PDE has smooth solutions (provided
$PV_{e}$-and-$M$ are smooth) which converge to the solution of the
unregularised $PV_{e}$-and-$M$ inversion. We will also observe that this
process is reversible: given an appropriately regularised version of
$PV_{e}$-and-$M$ inversion we may reconstruct the corresponding energy. We
then conclude this section with a discussion of the fact that minimizers of
the smoothed out energy converge to the minimizer of the true energy.

First we discuss the subtlety in determining the ellipticity of
$PV_{e}$-and-$M$ inversion.

######  Remark 6.1 (Agnostic smoothing fails to determine whether or not
$PV_{e}$-and-$M$ inversion is elliptic).

In order to streamline the computations we will begin with a purely algebraic
manipulation of $PV_{e}$-and-$M$ inversion in its form explicitly invoking
Heavisides. Since $H_{u}+H_{s}=1$, (1) may be rewritten as

$\Delta p+\frac{1}{2}\partial_{3}\left((M-\partial_{3}p)H_{u}\right)=PV_{e}.$
(21)

We now suppose that $\mathcal{H}_{u}$ is a smooth approximation of $H_{u}$,
remaining agnostic as to how precisely the approximation is made. We only
posit that $\mathcal{H}_{u}=\mathcal{H}_{u}(M,\,\partial_{3}p)$. This is
sensible since the original Heavisides depend on the water content $q$, which
itself may be written in terms of $M$ and $\partial_{3}p$ (the horizontal
derivatives of $p$ play no part in the determination of $q$). We thus consider

$\Delta
p+\frac{1}{2}\partial_{3}\left((M-\partial_{3}p)\mathcal{H}_{u}\right)=PV_{e}$
(22)

as a regularised version of (21). Since
$\mathcal{H}_{u}=\mathcal{H}_{u}(M,\,\partial_{3}p)$ is smooth we may compute
that

$\partial_{3}\left((M-\partial_{3}p)\mathcal{H}_{u}\right)=(\partial_{3}M-\partial_{3}^{2}p)\mathcal{H}_{u}+(M-\partial_{3}p)\left((\partial_{M}\mathcal{H}_{u})\partial_{3}M+(\partial_{p}\mathcal{H}_{u})\partial_{3}^{2}p\right)$

such that (22) may be written equivalently as

$\displaystyle
PV_{e}=\Delta_{h}p+\left(1-\frac{1}{2}\mathcal{H}_{u}+\frac{1}{2}(M-\partial_{3}p)\partial_{p}\mathcal{H}_{u}\right)\partial_{3}^{2}p$
$\displaystyle+\frac{1}{2}(\partial_{3}M)\mathcal{H}_{u}+\frac{1}{2}(M-\partial_{3}p)(\partial_{M}\mathcal{H}_{u})\partial_{3}M.$

The principal part is thus the operator

$\Delta_{h}+\left(1-\frac{1}{2}\mathcal{H}_{u}+\frac{1}{2}(M-\partial_{3}p)\partial_{p}\mathcal{H}_{u}\right)\partial_{3}^{2}.$

Here lies the difficulty: while sensible approximations would ensure that the
regularised Heaviside satisfies $0\leqslant\mathcal{H}_{u}\leqslant 1$ such
that $1-\frac{1}{2}\mathcal{H}_{u}$ is uniformly bounded below by $1/2$, we
know nothing about the sign of
$\frac{1}{2}(M-\partial_{3}p)\partial_{p}\mathcal{H}_{u}$, and therefore know
nothing about the sign of the coefficient of $\partial_{3}^{2}$. This means
that we cannot say anything about the ellipticity of $PV_{e}$-and-$M$
inversion with this approach.

We now turn our attention towards a more careful way to smooth out the
Heavisides, leveraging the variational formulation. This smoothing relies on
mollification – see Definition 4.2 for the convention used here – and at the
crux of the computations lies the expression for the smoothed out version of
the function $\min(\,\cdot\,,\,0)$ recorded below. This mollification is also
depicted in Figure 3.

######  Lemma 6.2 (Explicit formulae for ${\min}_{\varepsilon}$ and its
derivative).

We introduce, for any $x\in\mathbb{R}$, the notation
$\min_{0}(x)\vcentcolon=\min(x,\,0)$ for $x\in\mathbb{R}$ and, for
$\varphi_{\varepsilon}$ a standard mollifier,
$\min_{\varepsilon}\vcentcolon=\min_{0}\ast\varphi_{\varepsilon}$. If
$\varphi_{\varepsilon}$ is centered then

${\min}_{\varepsilon}(x)=\left\\{\begin{aligned} &x&&\text{if
}x\leqslant-\varepsilon,\\\ &F_{\varepsilon}(x)&&\text{if }\lvert
x\rvert<\varepsilon,\text{ and }\\\ &0&&\text{if
}x\geqslant\varepsilon\end{aligned}\right.$ (23)

where, for $\lvert x\rvert<\varepsilon$,

$F_{\varepsilon}(x)\vcentcolon=\frac{x}{2}+\frac{1}{2}\int_{x}^{\varepsilon}(x-y)\varphi_{\varepsilon}(y)dy-\frac{1}{2}\int_{-\varepsilon}^{x}(x-y)\varphi_{\varepsilon}(y)dy$
(24)

######  Proof 6.3.

We begin by noting that, since
$\operatorname*{supp}\varphi_{\varepsilon}\subseteq(-\varepsilon,\,\varepsilon)$,

${\min}_{\varepsilon}(x)=\int_{\max(-\varepsilon,\,x)}^{\varepsilon}(x-y)\varphi_{\varepsilon}(y)dy.$
(25)

To conclude the computation it thus suffices to split into three cases.

  * •

Case 1: $x\leqslant-\varepsilon$. Then $\max(-\varepsilon,\,x)=-\varepsilon$
and so, since $\varphi_{\varepsilon}$ is centered, (25) tells us that
$\min_{\varepsilon}(x)=x$.

  * •

Case 2: $x\geqslant\varepsilon$. Then
$\max(-\varepsilon,\,x)=x\geqslant\varepsilon$ and so the integral in (25) is
taken over an empty set, such that $\min_{\varepsilon}(x)=0$.

  * •

Case 3: $\lvert x\rvert<\varepsilon$ – a.k.a. the interesting case. Then
$\max(-\varepsilon,\,x)=x$. In particular since
$\int_{x}^{\varepsilon}=\int_{-\varepsilon}^{\varepsilon}-\int_{-\varepsilon}^{x}$
and since $\varphi_{\varepsilon}$ is centered we compute that

${\min}_{\varepsilon}(x)=\int_{x}^{\varepsilon}(x-y)\varphi_{\varepsilon}(y)dy=x-\int_{-\varepsilon}^{x}(x-y)\varphi_{\varepsilon}(y)dy.$

Symmetrizing these two expressions for $\min_{\varepsilon}$ produces, as
desired

${\min}_{\varepsilon}(x)=\frac{1}{2}\int_{x}^{\varepsilon}(x-y)\varphi_{\varepsilon}(y)dy+\frac{1}{2}\left(x-\int_{-\varepsilon}^{x}\varphi_{\varepsilon}(y)dy\right)=F_{\varepsilon}(x).$

The three cases above verify (23).

Figure 3:  Following the notation of Lemma 6.2, $\min_{\varepsilon}$ is the
mollification of $\min_{0}\vcentcolon=\min(\,\cdot\,,\,0)$. The explicit form
of the mollification is given in Lemma 6.2 in terms of $F_{\varepsilon}$ such
that $\min_{\varepsilon}=F_{\varepsilon}$ in $(-\varepsilon,\,\varepsilon)$
while $\min_{\varepsilon}=\min_{0}$ elsewhere. Here we show $\min_{0}$, the
reference mollification $F=F_{1}$, and $F_{\varepsilon}$ for
$0<\varepsilon<1$. Having an explicit expression for $\min_{\varepsilon}$ is
the key element required to have an explicit form of regularised
$PV_{e}$-and-$M$ inversion (see Theorem 6.4).

With the expressions for $\min_{\varepsilon}$ and its derivative obtained in
Lemma 6.2 above in hand we may deduce the form taken by the smoothed out
$PV_{e}$-and-$M$ inversion.

######  Theorem 6.4 (Regularised $PV_{e}$-and-$M$ inversion).

Let $\varphi_{\varepsilon}$ be a standard centered mollifier such that, for
$f_{0}(s)\vcentcolon=\frac{1}{2}{\min(s,\,0)}^{2}$ we define
$f_{\varepsilon}\vcentcolon=f_{0}\ast\varphi_{\varepsilon}$. Then for $M\in
L^{2}$, $PV_{e}\in H^{-1}$, and $p\in\mathring{H}^{1}$ we define

$E_{\varepsilon}(p)\vcentcolon=\int_{\mathbb{T}^{3}}\frac{1}{2}{\lvert\nabla
p\rvert}^{2}-\frac{1}{2}f_{\varepsilon}(M-\partial_{3}p)+\langle
PV_{e},\,p\rangle.$ (26)

For $\min_{\varepsilon}$ defined in Lemma 6.2 we have that, for any
$p,\,\phi\in\mathring{H}^{1}$,

$DE_{\varepsilon}(p)\phi=\int_{\mathbb{T}^{3}}\nabla
p\cdot\nabla\phi-\frac{1}{2}{\min}_{\varepsilon}(M-\partial_{3}p)(\partial_{3}\phi)+\langle
PV_{e},\,\phi\rangle$ (27)

such that $p\in\mathring{H}^{1}$ is the minimizer of $E_{\varepsilon}$ if and
only if $p$ is an $\mathring{H}^{1}$-weak solution of

$-\Delta
p-\frac{1}{2}\partial_{3}\left({\min}_{\varepsilon}(M-\partial_{3}p)\right)=-PV_{e}.$
(28)

In particular if we introduce

$H_{u}^{\varepsilon}\vcentcolon=\mathds{1}(M-\partial_{3}p\leqslant\varepsilon)\text{
and
}H_{int}^{\varepsilon}\vcentcolon=\mathds{1}(-\varepsilon<M-\partial_{3}p<\varepsilon)$
(29)

then we may write (28) as

$-\Delta
p-\frac{1}{2}\partial_{3}\left((M-\partial_{3}p)H_{u}^{\varepsilon}+F_{\varepsilon}(M-\partial_{3}p)H_{int}^{\varepsilon}\right),$
(30)

for $F_{\varepsilon}$ defined in (24).

######  Proof 6.5.

(27) is obtained from differentiating (26) since
$f_{\varepsilon}^{\prime}=f_{0}^{\prime}\ast\varphi_{\varepsilon}=\min_{\varepsilon}$.
Integrating by parts then yields (28) and using Lemma 6.2 produces yields
(30).

We now turn our attention to the matter of _reversibility_ : given an
appropriately smoothed out version of $PV_{e}$-and-$M$ inversion, can we
reconstruct the corresponding (smooth) energy? The answer is yes. In Lemma 6.6
we begin by recording some properties of the smoothed out version of
$PV_{e}$-and-$M$ inversion, then in Theorem 6.9 we show that if these
properties hold for _any_ smoothed out version of $PV_{e}$-and-$M$ inversion,
then a smooth energy can be reconstructed via mollification.

######  Lemma 6.6 (Properties of $F$).

Let $\varphi_{\varepsilon}$ be a standard centered mollifier and consider
$F_{\varepsilon}$ as in Lemma 6.2. Since $F_{\varepsilon}=\varepsilon
F\left(\frac{\cdot}{\varepsilon}\right)$ we define $F\vcentcolon=F_{1}$. Then
$F$ satisfies

$F(-1)=-1,\,F^{\prime}(-1)=1,\,F(1)=F^{\prime}(1)=0,\,F^{\prime\prime}\leqslant
0,\text{ and }\operatorname*{supp}F^{\prime\prime}\subseteq(-1,\,1).$ (31)

######  Proof 6.7.

The identities $F(-1)=-1$ and $F(1)=0$ follow immediately from $\varphi$ being
normalised and centered and from the definition of $F_{\varepsilon}$ in (24).
Similarly the identities $F^{\prime}(-1)=1$ and $F^{\prime}(1)=0$ follow from
$\varphi$ being normalised. Differentiating (24) twice tells us that
$F^{\prime\prime}=-\varphi$ and so that fact that $F^{\prime\prime}\leqslant
0$ follows from the non-negativity of $\varphi$.

######  Remark 6.8 (Demystifying the relation between $\varphi$ and $F$).

In the proof of Lemma 6.6 above we have obtained that
$F^{\prime\prime}=-\varphi$. This provides an alternative route to (24)
expressing $F$ in terms of $\varphi$. Indeed, by taking the mean of the
representations of $F(x)$ using a Taylor expansion about $-1$ and $1$, with
the remainder in integral form, tells us that

$\displaystyle
F(x)=\frac{F(-1)+F(1)}{2}+\frac{F^{\prime}(-1)-F^{\prime}(1)}{2}+\frac{F^{\prime}(-1)+F^{\prime}(1)}{2}x$
$\displaystyle+\frac{1}{2}\int_{-1}^{x}(x-y)F^{\prime\prime}(y)dy-\frac{1}{2}\int_{x}^{1}(x-y)F^{\prime\prime}(y)dy.$

Note that this identity holds for _any_ twice continuously differentiable
function. In particular, if (31) holds and $F^{\prime\prime}=-\varphi$ then we
recover (24).

We now state and prove the reversibility theorem.

######  Theorem 6.9 (Reversibility).

Suppose that a smooth function $F$ satisfies (31). Then there exists a
standard centered mollifier $\varphi$, given by $\varphi=-F^{\prime\prime}$,
such that if we define $E_{\varepsilon}$ as in (26), $H_{u}^{\varepsilon}$ and
$H_{\text{int}}^{\varepsilon}$ as in (29), and let
$F_{\varepsilon}\vcentcolon=\varepsilon
F\left(\frac{\cdot}{\varepsilon}\right)$ then, for any $p\in\mathring{H}^{1}$,
the strong form of $DE_{\varepsilon}(p)\phi=0$ for every
$\phi\in\mathring{H}^{1}$ is precisely (30).

######  Proof 6.10.

All that we need to do is verify that $\varphi$ is indeed a standard centered
mollifier. Clearly $\varphi=-F^{\prime\prime}$ is smooth, supported on
$(-1,\,1)$, and since $F^{\prime\prime}\leqslant 0$ by assumption we know that
$\varphi$ is non-negative. The fact that $\varphi$ is normalised follows from
the values of $F^{\prime}$ at $\pm 1$, and similarly the fact that $\varphi$
is centered follows from the identity $\int yF^{\prime\prime}=yF^{\prime}-\int
F^{\prime}$ and the values of $F$ and $F^{\prime}$ at $\pm 1$.

To conclude this section we turn our attention towards the proof that
minimizers of the smoothed out energy converge to the minimizer of the true
energy. In order to do so we first record the following result concerning the
derivative of the smoothed out energy.

######  Lemma 6.11 (The derivative of the $\varepsilon$-regularisation is
$\varepsilon$-close to the true derivative).

Let $M\in L^{2}$, $PV_{e}\in H^{-1}$, and let $E$ be the energy introduced in
Definition 4.1. Let $\varphi$ be a standard mollifier and let
$E_{\varepsilon}$ be as defined in Theorem 6.4. For any
$p\in\mathring{H}^{1}$, if we write $C_{\varphi}\vcentcolon=\int\lvert
y\rvert\varphi(y)dy$ then ${\lvert\lvert
DE(p)-DE_{\varepsilon}(p)\rvert\rvert}_{H^{-1}}\leqslant
C_{\varphi}{\lvert\mathbb{T}^{3}\rvert}^{1/2}\varepsilon/2.$

######  Proof 6.12.

As in Theorem 6.4 we write $f_{0}(s)\vcentcolon=\frac{1}{2}\min{(s,\,0)}^{2}$
and $f_{\varepsilon}\vcentcolon=f_{0}\ast\varepsilon$, noting that
$f_{0}^{\prime}=\min(\,\cdot\,,\,0)$. Subtracting the expression for
$DE_{\varepsilon}$ recorded in (28) from the expression for $DE$ recorded in
(7) then tells us that, for any $\phi\in\mathring{H}^{1}$,
$(DE(p)-DE_{\varepsilon}(p))\phi=\frac{1}{2}\int_{\mathbb{T}^{3}}(f_{0}^{\prime}-f_{\varepsilon}^{\prime})(M-\partial_{3}p)(\partial_{3}\phi).$
Since $\text{Lip}(f_{0}^{\prime})=1$ we then deduce from Hölder’s inequality
and Lemma B.25 that $\lvert(DE(p)-DE_{\varepsilon}(p))\phi\rvert\leqslant
C_{\varphi}\varepsilon{\lvert\mathbb{T}^{3}\rvert}^{1/2}{\lvert\lvert\phi\rvert\rvert}_{\mathring{H}^{1}}/2,$
from which the claim follows.

Finally we conclude this section by establishing that minimizers of the
smoothed out energy converge to the minimizer of the true energy.

######  Theorem 6.13 (Minimizers of the $\varepsilon$-energy are
$\varepsilon$-close to the minimizer of the true energy).

Let $M\in L^{2}$, $PV_{e}\in H^{-1}$, and let $E$ be the energy introduced in
Definition 4.1. Let $\varphi$ be a standard mollifier and let
$E_{\varepsilon}$ be as defined in Theorem 6.4. Denote
$p^{*}\vcentcolon=\operatorname*{arg\,min}_{\mathring{H}^{1}}E$ and
$p^{*}_{\varepsilon}=\operatorname*{arg\,min}_{\mathring{H}^{1}}E_{\varepsilon}$
and let $C_{\varphi}\vcentcolon=\int\lvert y\rvert\varphi(y)dy.$ Then
${\lvert\lvert
p^{*}-p^{*}_{\varepsilon}\rvert\rvert}_{\mathring{H}^{1}}\leqslant\sqrt{2}C_{\varphi}{\lvert\mathbb{T}^{3}\rvert}^{1/2}\varepsilon.$

######  Proof 6.14.

Since $DE(p^{*})=0$ the $\frac{1}{2}$–strong convexity of $E$, proved in
Proposition 1, in the form of (33) combined with Corollary 5 tells us that

${\lvert\lvert
p^{*}-p_{\varepsilon}^{*}\rvert\rvert}_{\mathring{H}^{1}}^{2}\leqslant
4\left[E(p^{*}\varepsilon)-E(p^{*})\right]\leqslant{\lvert\lvert
DE(P_{\varepsilon}^{*})\rvert\rvert}_{H^{-1}}^{2}.$

We may then use Lemma 6.11 to conclude since
$\nabla_{\mathring{H}^{1}}E_{\varepsilon}\left(p^{*}_{\varepsilon}\right)=0$
and hence, using the inequality above,

$\displaystyle{\lvert\lvert
p^{*}-p^{*}_{\varepsilon}\rvert\rvert}_{\mathring{H}^{1}}^{2}\leqslant
8\left({\lvert\lvert
DE_{\varepsilon}\left(p^{*}_{\varepsilon}\right)\rvert\rvert}_{H^{-1}}^{2}+{\lvert\lvert
DE_{\varepsilon}\left(p^{*}_{\varepsilon}\right)-DE\left(p^{*}_{\varepsilon}\right)\rvert\rvert}_{H^{-1}}^{2}\right)\leqslant
2C_{\varphi}^{2}\lvert\mathbb{T}^{3}\rvert\varepsilon^{2},$

as claimed.

##  Appendix A Tools from convex analysis

In this section we record various results from convex analysis that are of use
to us. This begins with results on various characterisations of strong
convexity and concludes with a lemma relating convexity to weak lower semi-
continuity of functionals.

First we recall a few first-order characterisations of strong convexity.

######  Lemma A.1 (Equivalent characterisations of strong convexity).

Let $\left(H,\,\langle\,\cdot\,,\,\cdot\,\rangle\right)$ be a Hilbert space
and let $f:H\to\mathbb{R}$ be Fréchet differentiable. For any $\mu>0$ the
following are equivalent, where each inequality holds for every $x,\,y\in H$
and $\theta\in[0,\,1]$.

$\displaystyle f(\theta x+(1-\theta)y)\leqslant\theta
f(x)+(1-\theta)f(y)-\theta(1-\theta)\frac{\mu}{2}{\lvert\lvert
x-y\rvert\rvert}^{2}$ (32) $\displaystyle f(y)\geqslant f(x)+\langle
Df(x),\,y-x\rangle+\frac{\mu}{2}{\lvert\lvert x-y\rvert\rvert}^{2}$ (33)
$\displaystyle\langle Df(x)-Df(y),\,x-y\rangle\geqslant\mu{\lvert\lvert
x-y\rvert\rvert}^{2}$ (34)

If any of these inequalities hold recall that we say that $f$ is _$\mu$
-strongly convex_.

######  Proof A.2.

See Definition 2.1.2 and Theorem 2.1.9 in [Nes04].

Similarly we now recall some second-order characterisations of strong
convexity.

######  Lemma A.3 (Second-order characterisation of strong convexity).

Let $\left(H,\,\langle\,\cdot\,,\,\cdot\,\rangle\right)$ be a Hilbert space
and let $f:H\to\mathbb{R}$ be twice continuously Fréchet differentiable. $f$
is $\mu$-strongly convex if and only if $D^{2}f\geqslant 0$, in the sense that
$D^{2}f(x)(\xi,\,\xi)\geqslant 0$ for every $x,\,\xi\in H$.

######  Proof A.4.

See Theorem 2.1.11 in [Nes04].

We now record a simple but very useful consequence of strong convexity which
tells us that the derivatives of strongly convex functions essentially act as
norms near their minimizers.

######  Corollary 5 (Derivative of a strongly convex function is essentially a
norm near its minimizer).

Let $\left(H,\,\langle\,\cdot\,,\,\cdot\,\rangle\right)$ be a Hilbert space
with dual $H^{*}$ and let $f:H\to\mathbb{R}$ be a Fréchet differentiable and
$\mu$-strongly convex. Then $f(x)-f(y)\leqslant\frac{1}{2\mu}{\lvert\lvert
Df(x)\rvert\rvert}_{H^{*}}^{2}$ for every $x,\,y\in H$.

######  Proof A.5.

We write $\nabla_{H}f(x)$ for the Riesz representative of $Df(x)$ and note
that $y^{*}=x-\frac{1}{2\mu}\nabla_{H}f(x)$ is the minimum of the right-hand
side of the characterisation (33) of strong convexity with respect to $y$.
Substituting $y=y^{*}$ into the right-hand side of (33) then tells us that,
for every $x,\,y\in H$,

$f(y)\geqslant
f(x)-\frac{1}{\mu}Df(x)\nabla_{H}f(x)+\frac{1}{2\mu}{\lvert\lvert\nabla_{H}f(x)\rvert\rvert}^{2}=f(x)-\frac{1}{2\mu}{\lvert\lvert
Df(x)\rvert\rvert}_{H^{*}}^{2}$

as desired since ${\lvert\lvert\nabla_{H}f(x)\rvert\rvert}={\lvert\lvert
Df(x)\rvert\rvert}_{H^{*}}$.

Below we now record a result used for the existence and uniqueness result of
Theorem 5.1, i.e. one of the ingredients of the Direct Methods of the Calculus
of Variations.

######  Lemma A.6 (Convex energy density implies weak lower semi-continuity).

Consider $f:\mathbb{T}^{3}\times\mathbb{R}^{3}\to\mathbb{R}$ which satisfies
the following conditions.

  1. 1.

For every $x\in\mathbb{T}^{3}$, $f(x,\,\cdot\,)$ is convex and continuously
differentiable.

  2. 2.

For every $u\in L^{2}(\mathbb{T}^{3},\,\mathbb{R}^{3})$,

$x\mapsto f(x,\,u(x))\in L^{1}(\mathbb{T}^{3})\text{ and
}x\mapsto(\nabla_{u}f)(x,\,u(x))\in L^{2}(\mathbb{T}^{3}).$

The functional $\mathcal{F}:\mathring{H}^{1}\to\mathbb{R}$ defined by
$\mathcal{F}(p)\vcentcolon=\int_{\mathbb{T}^{3}}f(x,\,\nabla p(x))dx$ for
every $p\in\mathring{H}^{1}$ is weakly lower semi-continuous in
$\mathring{H}^{1}$.

######  Proof A.7.

Let $(p_{n})$ be a sequence in $\mathring{H}^{1}$ weakly convergent to
$p\in\mathring{H}^{1}$. Since $f(x,\,\cdot\,)$ is convex and continuously
differentiable we have that

$f(x,\,\nabla p)\leqslant f(x,\,\nabla p_{n})+(\nabla_{u}f)(x,\,\nabla
p)\cdot(\nabla p-\nabla p_{n}).$

In particular, since $(\nabla_{u}f)(x,\,\nabla p)\in L^{2}$ and $\nabla
p-\nabla p_{n}\rightharpoonup 0\text{ in }L^{2},$ we deduce that indeed

$\mathcal{F}(p)=\int_{\mathbb{T}^{3}}f(x,\,\nabla
p)\leqslant\liminf\int_{\mathbb{T}^{3}}f(x,\,\nabla
p_{n})+0=\liminf\mathcal{F}(p_{n}),$

as desired.

##  Appendix B Tools from regularity theory

In this section we record tools from regularity theory that are of use to us.
First we discuss elementary properties of finite differences (used to obtain
the $H^{2}$ regularity of solutions in Lemma 5.3). The bulk of this section is
then devoted to the discussion of the de Giorgi regularity theory, and
specifically its formulation in terms of Campanato functions. Finally this
section concludes with a simple result pertaining to approximations using
mollifiers.

First we fix notation and recall the definition of finite differences.

######  Definition B.1 (Finite differences).

We define, for $f\in L^{1}_{\text{loc}}(\mathbb{T}^{3})$,
$x,\,h\in\mathbb{T}^{3}$, and $i=1,\,\dots,\,n$ the _finite difference_
$\Delta_{h}^{i}f(x)\vcentcolon=f(x+he_{i})-f(x).$

We now record elementary properties of finite differences.

######  Lemma B.2 (Properties of finite differences).

For finite differences the following hold.

  1. 1.

Adjoint:
$\int_{\mathbb{T}^{3}}(\Delta_{h}^{i}f)g=\int_{\mathbb{T}^{3}}f(\Delta_{-h}^{i}g)$.

  2. 2.

Commutation: $\partial_{j}$ and $\Delta_{h}^{i}$ commute.

  3. 3.

Chain Rule: if $f$ is absolutely continuous then the following identity holds:
$\Delta_{h}^{i}(f\circ
g)=(\Delta_{h}^{i}g)\int_{0}^{1}f^{\prime}\left((\operatorname{id}+\theta\Delta_{h}^{i})g\right)d\theta.$

  4. 4.

Characterisation of $W^{1,\,p}$, necessity: If $f\in
W^{1,\,p}(\mathbb{T}^{3})$ for $p\in[1,\,\infty)$ then we have that
${\lvert\lvert\Delta_{h}^{i}f\rvert\rvert}_{L^{p}}\leqslant C_{d}\lvert
h\rvert\,{\lvert\lvert\nabla f\rvert\rvert}_{L^{p}}.$

  5. 5.

Characterisation of $W^{1,\,p}$, sufficiency: If $u\in L^{p}(\mathbb{T}^{3})$,
$p\in(1,\,\infty)$, such that
${\lvert\lvert\Delta_{h}^{i}u\rvert\rvert}_{L^{p}}\leqslant M\lvert h\rvert$
for some constant $M>0$ then $u\in W^{1,\,p}(\mathbb{T}^{3})$, we have the
estimate ${\lvert\lvert\nabla u\rvert\rvert}_{L^{p}}\leqslant M$, and
$\frac{1}{h}\Delta_{h}^{i}u\rightharpoonup\partial_{i}u$ in $L^{p}$ as $h\to
0$.

  6. 6.

$H^{-1}$ estimate: If $u\in L^{2}$ then $\Delta_{h}^{i}f\in H^{-1}$ and
${\lvert\lvert\Delta_{h}^{i}f\rvert\rvert}_{H^{-1}}\leqslant\lvert
h\rvert\,{\lvert\lvert f\rvert\rvert}_{L^{2}}.$

######  Proof B.3.

Item 1 follows from a direct computation. Item 2 follows from the fact that
$\Delta_{h}^{i}f=f\circ\Phi_{h}^{i}-f$ for
$\Phi_{h}^{i}(x)\vcentcolon=x+he_{i}$, such that $D\Phi_{h}^{i}=I$. Items 4
and 5 are proved in [GM12]. Item 6 follows from the fact that, for any $f\in
L^{2}$, $\nabla f\in H^{-1}$ with ${\lvert\lvert\nabla
f\rvert\rvert}_{H^{-1}}\leqslant{\lvert\lvert f\rvert\rvert}_{L^{2}}$. Now we
finally we turn our attention to the proof of item 3, which comes down to a
short computation. Writing $y\vcentcolon=g(x)$ and $z\vcentcolon=g(x+he_{i})$,
such that $z-y=(\Delta_{h}^{i}g)(x)$, we have that

$\displaystyle\left[\Delta_{h}^{i}(f\circ
g)\right](x)=f(z)-f(y)=\int_{0}^{1}(z-y)f^{\prime}\left(y+\theta(z-y)\right)d\theta$
$\displaystyle=(\Delta_{h}^{i}g)(x)\int_{0}^{1}f^{\prime}\left(g(x)+\theta(\Delta_{h}^{i}g)(x)\right)d\theta,$

as claimed.

We now turn our attention towards the topic occupying the bulk of this
section: formulating the de Giorgi regularity theory using Campanato
functions. There are three important components of that in this section: (1)
obtaining a sufficient condition for Hölder continuity in terms of local
integrals of gradients, which is done in Proposition 6, (2) reformulating the
classical de Giorgi result in terms of Campanato functions, which is done in
Proposition 8, and (3) deriving an _a priori_ estimate of Campanato-type for
elliptic equations with forcing in divergence-form.

First we fix notation and define Campanato functions.

######  Definition B.4 (Campanato functions).

We say that $u\in L^{2}(\mathbb{T}^{d})$ is _Campanato (regular) of order
$\alpha$_ for some $\alpha\in(0,\,1)$ if there exists a constant $K>0$ such
that, for every $x\in\mathbb{R}^{3}$ and $R>0$, $\int_{B(x,\,R)}{\lvert
u-u_{x,\,R}\rvert}^{2}\leqslant K^{2}R^{d+2\alpha}$ where
$u_{x,\,R}\vcentcolon=\fint_{B(x,\,R)}u$ denotes the local average of $u$. We
call $K$ the _Campanato constant_ of $u$.

######  Remark B.5 (Interpretation of the Campanato condition).

It is interesting to consider the form of the defining estimate of Campanato
functions using _normalised averages_ and we do this here in the (slightly
more general) context of $L^{p}$ based Campanato spaces, for
$p\in[1,\,\infty)$. Specifically, note that the defining inequality
$\int_{B(x,\,R)}{\lvert u-u_{x,\,R}\rvert}^{2}\leqslant M^{p}R^{d+\alpha p}$
is equivalent to ${\left(\fint_{B(x,\,R)}{\lvert
u-\overline{u}_{x,\,R}\rvert}^{p}\right)}^{1/p}\leqslant\frac{M}{\omega_{d}^{1/p}}R^{\alpha},$
where $\omega_{d}$ denotes the $d$-dimensional volume of the unit ball in
$\mathbb{R}^{d}$. The latter inequality makes it more apparent that the
deviation of $u$ from its average is only controlled by $R^{\alpha}$, as would
be expected from a function which is $\alpha$-Hölder continuous.

We now record half of the key fact about Campanato function, which is that
they are nothing more than an equivalent characterisation of Hölder continuous
functions. We only mention one direction of the equivalence below, the other
direction is much easier to prove.

######  Lemma B.6 (Campanato implies Hölder).

Let $u\in L^{2}(\mathbb{T}^{d})$ be Campanato of order $\alpha$, for some
$\alpha\in(0,\,1)$, with Campanato constant $K$. Then $u$ is Hölder continuous
with exponent $\alpha$, and for every open set
$\mathcal{U}\subseteq\mathbb{T}^{3}$, ${\lvert\lvert
u\rvert\rvert}_{C^{\alpha}(\mathcal{U})}\leqslant C\left(K+{\lvert\lvert
u\rvert\rvert}_{L^{2}}\right)$ for some constant
$C=C(\alpha,\,d,\,\mathcal{U})$.

######  Proof B.7.

See Theorem 3.1 of [HL11]

We now progress towards the aforementioned sufficient condition for Hölder
continuity involving local integrals of gradients. The missing ingredient is a
careful analysis of the scaling of the Poincaré constant in balls.

######  Lemma B.8 (Scaling of the Poincaré constant in balls).

There exists a constant $C=C(d)>0$ such that, for every $x\in\mathbb{T}^{d}$,
every $r>0$, and every $u\in H^{1}((B(x,\,r))$, the inequality
$\int_{B(x,\,r)}{\big{|}u-\fint_{B(x,\,r)}u\big{\rvert}}^{2}\leqslant
Cr^{2}\int_{B(x,\,r)}{\lvert\nabla u\rvert}^{2}$ holds.

######  Proof B.9.

This follows from the Poincaré inequality in the unit ball and scaling:
defining $v(y)\vcentcolon=u(x+ry)$ such that $v\in H^{1}(B(0,1))$ and applying
the Poincaré inequality to $v$ produces the claim.

We are now ready to state and prove the first of three main results of this
section.

######  Proposition 6 (Sufficient Campanato condition involving local
integrals of gradients).

Suppose that $u\in H^{1}_{\text{loc}}(\mathbb{T}^{d})$ such that, for some
$K>0$ and $\alpha\in(0,\,1)$, the inequality $\int_{B(x,\,r)}{\lvert\nabla
u\rvert}^{2}\leqslant K^{2}r^{d+2\alpha-2}$ holds for every
$B(x,\,r)\subseteq\mathbb{T}^{d}$. Then $u\in C^{\alpha}$ and, for every open
set $\mathcal{U}\subseteq\mathbb{T}^{d}$, we have the estimate ${\lvert\lvert
u\rvert\rvert}_{C^{\alpha}(\mathcal{U})}\leqslant C\left(K+{\lvert\lvert
u\rvert\rvert}_{L^{2}(\mathbb{T}^{d})}\right)$ for some constant
$C=C\left(\alpha,\,d,\,\mathcal{U}\right)$.

######  Proof B.10.

This follows from combining Lemma B.6, which says that all Campanato functions
are Hölder, with the Poincaré inequality in a ball – c.f. Lemma B.8.

We now turn our attention towards proving the second of three main results of
this section, namely reformulating the classical de Giorgi Theorem in terms of
Campanato functions. First we need this technical lemma about test functions,
which really tells us that we can find a smooth version $\varphi$ of the
indicator $\mathds{1}_{B_{r}}$ which vanishes outside the larger ball $B_{R}$
while retaining quantitative control over the gradient of $\varphi$.

######  Lemma B.11 (Careful scaling of a test function).

There exists a universal constant $C>0$ such that for every $0<r<R$ we may
find $\varphi\in C_{c}^{\infty}(B_{R})$ satisfying $0\leqslant\varphi\leqslant
1$ such that $\varphi\equiv 1$ in $B_{r}$ and
${\lvert\lvert\nabla\varphi\rvert\rvert}_{\infty}\leqslant\frac{C}{R-r}$.

######  Proof B.12.

Consider a radial profile $\rho\in C^{\infty}(\mathbb{R})$ satisfying
$0\leqslant\rho\leqslant 1$ which is supported in $(-\infty,\,1)$ and such
that $\rho(s)=1$ if $s<0$. The function
$\varphi(x)\vcentcolon=\rho\left(\frac{\lvert x\rvert-r}{R-r}\right)$
satisfies the desired properties, with
$C={\lvert\lvert\rho^{\prime}\rvert\rvert}_{\infty}$.

Since we will refer to uniformly elliptic matrices repeatedly in the sequel we
fix notation on the matter here.

######  Definition B.13 (Uniform ellipticity).

Let $\mathcal{U}\subseteq\mathbb{T}^{3}$. A matrix-valued pointwise symmetric
function $A\in L^{\infty}$ is said to be _uniformly elliptic with constants
$\theta,\,\Lambda>0$ in $\mathcal{U}$_ if, for every $x\in\mathcal{U}$,
$A(x)\geqslant\theta I$ and $\lvert A(x)\rvert\leqslant\Lambda$ for every
$x\in\mathbb{T}^{3}$.

We now continue building towards the Campanato formulation of the classical de
Giorgi result – the next step is to obtain the Caccioppoli inequality. It will
follow from the following energy estimate.

######  Lemma B.14 (Energy inequality, a.k.a. prelude to the Caccioppoli
inequality).

Let $u\in H^{1}$ be an $H^{1}$-weak solution of $-\nabla\cdot\left(A(x)\nabla
u\right)=0$ in $\mathbb{T}^{d}$ where $A\in L^{\infty}$ is pointwise symmetric
and uniformly elliptic with constants $\theta,\,\Lambda>0$ in
$\mathbb{T}^{d}$. For every $\varphi\in C_{c}^{\infty}$ the following estimate
holds: $\int\varphi^{2}{\lvert\nabla
u\rvert}^{2}\leqslant\frac{4\Lambda^{2}}{\theta^{2}}\int
u^{2}{\lvert\nabla\varphi\rvert}^{2}.$

######  Proof B.15.

We may test $-\nabla\cdot(A\nabla u)=0$ against $\varphi^{2}u$ to observe that
the identity $\int_{\mathbb{T}^{3}}\varphi^{2}A\nabla u\cdot\nabla
u=-2\int_{\mathbb{T}^{3}}A(\varphi\nabla u)\cdot(u\nabla\varphi)$ holds. Using
the ellipticity of $A$ on the left-hand side and an $\varepsilon$-Cauchy
inequality on the right-hand side of this identity then allows us to conclude.

We are now equipped to prove the Caccioppoli inequality.

######  Corollary 7 (Caccioppoli inequality).

Let $u\in H^{1}$ be an $H^{1}$-weak solution of $-\nabla\cdot\left(A(x)\nabla
u\right)=0$ in $\mathbb{T}^{d}$ where $A\in L^{\infty}$ is pointwise symmetric
and uniformly elliptic with constants $\theta,\,\Lambda>0$ in
$\mathbb{T}^{d}$. There exists a constant $C=C(d)>0$ such that, for every
$0<r<R<\infty$ and every $a\in\mathbb{R}$, $\int_{B_{r}}{\lvert\nabla
u\rvert}^{2}\leqslant\frac{C\Lambda^{2}}{{\theta^{2}\left(R-r\right)}^{2}}\int_{B_{R}}{\lvert
u-a\rvert}^{2}.$

######  Proof B.16.

Without loss of generality we may set $m=0$ since otherwise $v\vcentcolon=u-a$
satisfies $\nabla v=\nabla u$ and hence $v$ solves the same equation as $u$.
We then pick $\varphi$ as in Lemma B.11 and use it as a test function in Lemma
B.14 to establish the claim.

We are now ready to reformulate the classical de Giorgi result in the language
of Campanato. First we record that classical theorem here.

######  Theorem B.17 (de Giorgi).

Let $u\in H^{1}(\Omega)$ be an $H^{1}$-weak solution of the equation
$-\nabla\cdot\left(A(x)\nabla u\right)=0$ in a domain
$\Omega\subseteq\mathbb{R}^{d}$ with $d\geqslant 3$, where $A\in L^{\infty}$
is pointwise symmetric and uniformly elliptic with constants
$\theta,\,\Lambda>0$ in $\Omega$. There exists
$\alpha=\alpha(d,\,\theta,\,\Lambda)\in(0,\,1)$ and
$C=C(d,\,\theta,\,\Lambda)>0$ such that $u\in C^{\alpha}$ and, for every
$x_{0}$,

$\sup_{x\in B\left(x_{0},\,1/4\right)}\frac{\lvert
u(x)-u(x_{0})\rvert}{{\lvert x-x_{0}\rvert}^{\alpha}}\leqslant C{\lvert\lvert
u\rvert\rvert}_{L^{2}(B(x_{0},\,2))}.$

######  Proof B.18.

This was first proved in [DG57]. Very digestible proofs are provided in [CV10,
Vas16].

When reformulated using Campanato functions, de Giorgi’s Theorem becomes the
following result.

######  Proposition 8 (Campanato formulation of de Giorgi’s Theorem).

Let the function $w\in H^{1}(B(x_{0},\,R))\subseteq\mathbb{T}^{d}$ be an
$H^{1}$-weak solution of $-\nabla\cdot\left(A(x)\nabla w\right)=0$ in
$B(x_{0},\,R)$ where $A\in L^{\infty}$ is pointwise symmetric and uniformly
elliptic with constants $\theta,\,\Lambda$ in $B(x_{0},\,R)$. There exist
$\alpha=\alpha(d,\,\theta,\,\Lambda)\in(0,\,1)$ and
$C=C(d,\,\theta,\,\Lambda)>0$ such that, for any $0<\rho<R$,
$\int_{B(x_{0},\,\rho)}{\lvert\nabla w\rvert}^{2}\leqslant
C{\left(\frac{\rho}{R}\right)}^{d-2+2\alpha}\int_{B(x_{0},\,R)}{\lvert\nabla
w\rvert}^{2}.$

######  Proof B.19.

Note that if the claim is proved for $x_{0}=0$ and $R=2$ then we may deduce
that it holds for arbitrary $x_{0}$ and $R$ by translation and scaling.
Indeed, we may define $u(x)\vcentcolon=w\left(x_{0}+\frac{R}{2}x\right)$ such
that $u$ solves $-\nabla\cdot(\widetilde{A}\nabla u)=0$ in $B(0,\,2)$, where
crucially $\widetilde{A}(x)\vcentcolon=A\left(x_{0}+\frac{R}{2}x\right)$ has
the same uniform ellipticity constants as $A$. This allows us to deduce that
if the claim holds for $w$ in $B(0,\,2)$, then it holds for $u$ in
$B(x_{0},\,R)$.

So now we show that the claim holds when $x_{0}=0$ and $R=2$. Since $u-\fint
u$ solves the same equation as $u$ we may without loss of generality assume
that $\fint_{B_{2}}u=0$. We now combine Corollary 7, Theorem B.17, and the
Poincaré inequality:

$\displaystyle\int_{B_{\rho}}{\lvert\nabla
u\rvert}^{2}\leqslant\frac{C}{\rho^{2}}\int_{B_{2\rho}}{\lvert
u-u(0)\rvert}^{2}\leqslant\frac{C}{\rho^{2}}\int_{B_{2\rho}}{\lvert
x\rvert}^{2\alpha}\left(\int_{B_{2}}u^{2}\right)\leqslant
C\rho^{d+2\alpha-2}\int_{B_{2}}{\lvert\nabla u\rvert}^{2},$

where we have used that $\int_{B_{r}}{\lvert x\rvert}^{p}=C(d,\,p)s^{d+p}$.
Note that the chain of inequalities above only holds if $\rho\leqslant 1/8$,
such that $2\rho\leqslant 1/4$ and Theorem B.17 applies. Nonetheless, the
desired inequality also holds for $\rho>1/8$, albeit more trivially so since
in that case we only need to note that $\rho^{d+2\alpha-2}\geqslant
8^{2-2\alpha-d}=C(\alpha,\,d)$, and this concludes the proof.

We now record the third of three main results from this section: an _a priori_
estimate for elliptic equations whose forcing is in divergence-form.

######  Proposition 9 (Campanato-type _a priori_ estimate for forcing in
divergence form).

Suppose that $d>1$ and let $w\in H_{0}^{1}(B(x_{0},\,R))$ be an
$H_{0}^{1}$-weak solution of $-\nabla\cdot\left(A\nabla w\right)=\nabla\cdot
F$ in $B(x_{0},\,R)\subseteq\mathbb{T}^{d}$, where $A\in L^{\infty}$ is
pointwise symmetric and uniformly elliptic with constants $\theta,\,\Lambda>0$
in $\mathbb{T}^{d}$, and where $F\in L^{q}(B(x_{0},\,R))$ for some $q>d$.
There exist $\alpha=\alpha(d,\,q)\in(0,\,1)$ and $C=C(\alpha,\,d,\,\theta)>0$
such that $\int_{B(x_{0},\,R)}{\lvert\nabla w\rvert}^{2}\leqslant
C{\lvert\lvert F\rvert\rvert}_{L^{q}(B(x_{0},\,R))}R^{d-2+2\alpha}.$

######  Proof B.20.

We choose $\tilde{q}\in(1,\,\infty)$ such that
$\frac{1}{q}+\frac{1}{\tilde{q}}=\frac{1}{2}$ , which is possible since
$q>d\geqslant 2$. The estimate then follows from combining the ellipticity of
$A$, integration by parts, and the Hölder inequality (twice) since, using that
$w=0$ on $\partial B(x_{0},\,R)$,

$\displaystyle\theta{\lvert\lvert\nabla
w\rvert\rvert}_{L^{2}(B(x_{0},\,R))}^{2}$
$\displaystyle\leqslant\int_{B(x_{0},\,R)}A\nabla w\cdot\nabla
w=-\int_{B(x_{0},\,R)}F\cdot\nabla w$ $\displaystyle\leqslant{\lvert\lvert
F\rvert\rvert}_{L^{q}(B(x_{0},\,R))}{\lvert
B(x_{0},\,R)\rvert}^{d/\tilde{q}}{\lvert\lvert\nabla
w\rvert\rvert}_{L^{2}(B(x_{0},\,R))}.$

First dividing both sides of this inequality by ${\lvert\lvert\nabla
w\rvert\rvert}_{L^{2}}^{2}$ and then squaring yields ${\lvert\lvert\nabla
w\rvert\rvert}_{L^{2}(B(x_{0},\,R)}^{2}\leqslant
C(d,\,q,\,\theta){\lvert\lvert
F\rvert\rvert}_{L^{q}(B(x_{0},\,R)}^{2}R^{2d/\tilde{q}}.$ To conclude we
simply note that, since $q>d$, $\alpha\vcentcolon=1-\frac{d}{q}$ belongs to
$(0,\,1)$ such that $2d/\tilde{q}=d-2+2\alpha$ as desired.

As we approach the close of the portion of this section devoted to results
pertaining to the de Giorgi theory we record the following technical lemma
without proof.

######  Lemma B.21.

Let $\phi:(0,\,\infty)\to(0,\,\infty)$ be a non-decreasing map. Suppose that
there are $0<\beta<\gamma$ and $A,\,B,\,R_{0}>0$ such that
$\phi(\rho)\leqslant
A{\left(\frac{\rho}{r}\right)}^{\gamma}\phi(r)+Br^{\beta}$ for every
$0<\rho\leqslant r\leqslant R_{0}$. Then, for every
$\delta\in(\beta,\,\gamma)$ there exists $C=C(A,\,\beta,\,\gamma,\,\delta)>0$
such that $\phi(\rho)\leqslant
C\left[{\left(\frac{\rho}{r}\right)}^{\delta}\phi(r)+B\rho^{\beta}\right]$ for
every $0<\rho\leqslant r\leqslant R_{0}$.

######  Proof B.22.

See Lemma 3.4 of [HL11].

The last result in this section relevant to the de Giorgi theory is the
following, which is a standard elliptic estimate.

######  Lemma B.23 (Standard $\mathring{H}^{1}$ estimate for uniformly
elliptic operators).

Let $u$ be an $\mathring{H}^{1}$-weak solution of $-\nabla\cdot(A(x)\nabla
u)=\nabla\cdot F$ in $\mathbb{T}^{d}$ for some $A\in L^{\infty}$ which is
uniformly elliptic with constants $\theta,\,\Lambda$ in $\mathbb{T}^{d}$ and
some $F\in L^{2}(\mathbb{T}^{d})$. The following estimate holds:
${\lvert\lvert\nabla
u\rvert\rvert}_{L^{2}}\leqslant\frac{2}{\theta}{\lvert\lvert
F\rvert\rvert}_{L^{2}}$.

######  Proof B.24.

This follows from a standard energy estimate, using $u$ as a test function:

$\displaystyle\theta\int_{\mathbb{T}^{d}}{\lvert\nabla
u\rvert}^{2}\leqslant\int_{\mathbb{T}^{d}}A(x)\nabla u\cdot\nabla
u=\int_{\mathbb{T}^{d}}F(\nabla\cdot u)\leqslant\frac{2}{\theta}{\lvert\lvert
F\rvert\rvert}_{L^{2}}^{2}+\frac{\theta}{2}{\lvert\lvert\nabla
u\rvert\rvert}_{L^{2}}^{2}.$

An absorbing $\varepsilon$-Cauchy inequality concludes the proof.

Finally we conclude this section with a convergence result for approximations
using mollifiers which is used in Section 6.

######  Lemma B.25 (Uniform convergence of mollification of Lipschitz
functions).

Let $\varphi$ be a standard mollifier and let $f:\mathbb{R}\to\mathbb{R}$ be
Lipschitz. For $C_{\varphi}\vcentcolon=\int\lvert y\rvert\varphi(y)dy<\infty,$
we have that
$\lvert\left(f-f\ast\varphi_{\varepsilon}\right)(x)\rvert\leqslant
C_{\varphi}\text{Lip}(f)\varepsilon$ for any $x\in\mathbb{R}$ and all
$\varepsilon>0$.

######  Proof B.26.

Note that
$f(x)-(f\ast\varphi_{\varepsilon})(x)=\int\left[f(x)-f(x-y)\right]\varphi_{\varepsilon}(y)dy$
since $\varphi$ is normalised. The claim then follows from the non-negativity
of $\varphi$ and the change of variable $z=y/\varepsilon$.

##  Appendix C Rosetta stone

The purpose of this section is to provide a dictionary allowing us to
translate between the notation of [SS17] which introduces the precipitating
quasi-geostrophic system and the notation of this paper. In particular we will
discuss how considering the simple version of $PV_{e}$-and-$M$ inversion
recorded in (1) may be done without loss of generality, qualitatively
speaking.

For the purpose of this section, and of this section only, we introduce the
constants $f$, $B_{\theta}$, $B_{q}$, $C_{\theta}$, and $C_{q}>0$. We warn the
reader that they should not get too attached to these constants: they only
appear in this section and are set to unity (and thus disappear) everywhere
else in the paper. They appear here for one reason: justifying that the
version of $PV_{e}$-and-$M$ inversion recorded in (1) is not _overly_
simplified.

Table 1 provides a dictionary between the constants introduced above and those
used in [SS17], as well as between the notation used here and the notation
used in that paper.

| 2017 paper | Our notation
---|---
$\hat{z}$, $\partial_{z}$ |  $e_{3}$, $\partial_{3}$
$\nabla_{h}$ | $\nabla_{h}$
$(u,\,v,\,0)$ | $u$
$\zeta$ | ${(\nabla\times u)}_{3}$
$\phi$ | $p$
$\psi$ | $p/f$
$\theta_{e}$ | $\theta_{e}$
$q_{t}$ | $q$
$PV_{e}$ | $PV_{e}$
$M$ | $M$
| 2017 paper | Our notation
---|---
$f$ | $f$
$\frac{g}{\theta_{0}}$ | $B_{\theta_{e}}$
$\frac{gL_{v}}{\theta_{0}c_{p}}$ | $B_{q}$
$\frac{d\tilde{\theta}_{e}}{dz}$ | $C_{\theta_{e}}$
$-\frac{d\tilde{q}_{t}}{dz}=-\frac{d\tilde{q}_{v}}{dz}$ | $C_{q}$
(as $\tilde{q}_{t}=\tilde{q}_{v}$) |
Table 1:  Each table compares the notation used in [SS17] (in its left column)
and the notation used in this paper (in its right column). Note that the
constants appearing in the right column of the table on the right-hand side
are only used in Appendix C and are all set to unity elsewhere.

With this dictionary in hand we can compare the forms taken by various
expressions and identities, depending on the notation used. This is recorded
in Table 2, which contains expressions for geostrophic and hydrostatic
balance, as well as for the potential vorticity, the moist variable $M$, and
the saturated and unsaturated buoyancy frequencies (introduced in [SS17]).

2017 paper | Our notation | Constants set to unity
---|---|---
$f\hat{z}\times(u,\,v,\,0)=-\nabla_{h}\phi$ | $u=\frac{1}{f}\nabla_{h}^{\perp}h$ | $u=\nabla_{h}^{\perp}p$
$f\frac{\partial\psi}{\partial z}=\frac{g}{\theta_{0}}\theta_{e}-\frac{gL_{v}}{\theta_{0}c_{p}}q_{t}H_{u}$ | $\partial_{3}p=B_{\theta_{e}}\theta_{e}-B_{q}qH_{u}$ | $\partial_{3}p=\theta_{e}-qH_{u}$
$PV_{e}=\zeta+\frac{f}{d\tilde{\theta}_{e}/dz}\frac{\partial\theta_{e}}{\partial z}$ | $PV_{e}={(\nabla\times u)}_{3}+\frac{f}{C_{\theta_{e}}}\partial_{3}\theta_{e}$ | $PV_{e}={(\nabla\times u)}_{3}+\partial_{3}\theta_{e}$
$M=q_{t}-\frac{d\tilde{q}_{t}/dz}{d\tilde{\theta}_{e}/dz}\theta_{e}$ | $M=q+\frac{C_{q}}{C_{\theta_{e}}}\theta_{e}$ | $M=\theta_{e}+q$
$N_{u}^{2}=\frac{g}{\theta_{0}}\frac{d\tilde{\theta}_{e}}{dz}-\frac{gL_{v}}{\theta_{0}c_{p}}\frac{d\tilde{q}_{t}}{dz}$ | $N_{u}^{2}=B_{\theta_{e}}C_{\theta_{e}}+B_{q}C_{q}$ | $N_{u}^{2}=2$
$N_{s}^{2}=\frac{g}{\theta_{0}}\frac{d\tilde{\theta}_{e}}{dz}$ | $N_{s}^{2}=B_{\theta_{e}}C_{\theta_{e}}$ | $N_{s}^{2}=1$
Table 2:  Comparisons of various identities and expressions using the notation
of [SS17] (on the left), the notation of Appendix C (in the middle), and the
notation of the everywhere else in this paper (on the right). From top to
bottom: geostrophic balance, hydrostatic balance, (equivalent) potential
vorticity, moist variable $M$, the unsaturated buoyancy frequency, and the
saturated buoyancy frequency.

With all of these notational equivalences in hand, we may formulate
$PV_{e}$-and-$M$ inversion in six different ways: two for each of the three
sets of notation; one with the Heavisides $H_{u}$ and $H_{s}$ and one invoking
the function $\min_{0}=\min(\,\cdot\,,\,0)$. From top to bottom we use the
notation of [SS17], the notation of this section, and the notation used
everywhere else in this paper (which is the same as the notation used in this
section, albeit with the constants $f$, $B_{\theta_{e}}$, $B_{q}$,
$C_{\theta_{e}}$, and $C_{q}$ set to unity):

$\displaystyle PV_{e}$
$\displaystyle=\nabla^{2}_{h}\psi+f^{2}\frac{\partial}{\partial
z}\left[\frac{1}{N_{u}^{2}}\left(\frac{\partial\psi}{\partial
z}+\frac{gL_{v}}{\theta_{0}fc_{p}}M\right)H_{u}+\frac{1}{N_{s}^{2}}\left(\frac{\partial\psi}{\partial
z}\right)H_{s}\right]$
$\displaystyle=\nabla^{2}_{h}\psi+f^{2}\frac{\partial}{\partial
z}\left[\frac{1}{N_{s}^{2}}\frac{\partial\psi}{\partial
z}+\left(\frac{1}{N_{s}^{2}}-\frac{1}{N_{u}^{2}}\right){\min}_{0}\,\left(\frac{-g}{f\theta_{0}}\cdot\frac{d\tilde{\theta}_{e}/dz}{d\tilde{q}_{t}/dz}M-\partial_{3}p\right)\right]$
$\displaystyle PV_{e}$
$\displaystyle=\frac{1}{f}\Delta_{h}p+f\left[\frac{1}{N_{u}^{2}}\left(\partial_{3}p+B_{q}M\right)H_{u}+\frac{1}{N_{s}^{2}}\left(\partial_{3}p\right)H_{s}\right]$
(35)
$\displaystyle=\frac{1}{f}\Delta_{h}p+f\partial_{3}\left[\frac{1}{N_{s}^{2}}\partial_{3}p+\left(\frac{1}{N_{s}^{2}}-\frac{1}{N_{u}^{2}}\right){\min}_{0}\,\left(\frac{B_{\theta_{e}}C_{\theta_{e}}}{C_{q}}M-\partial_{3}p\right)\right]$
$\displaystyle PV_{e}$
$\displaystyle=\Delta_{h}p+\partial_{3}\left[\frac{1}{2}(\partial_{3}p+M)H_{u}+(\partial_{3}p)H_{s}\right]$
$\displaystyle=\Delta_{h}p+\partial_{3}\left[\partial_{3}p+\frac{1}{2}{\min}_{0}\,\left(M-\partial_{3}p\right)\right]=\Delta
p+\frac{1}{2}\partial_{3}{\min}_{0}\,\left(M-\partial_{3}p\right)$

In particular we note that the equations above follow from the dictionary
recorded in Table 2. The only intermediate computation required is recorded
below, without proof.

######  Lemma C.1.

Suppose that, for some strictly positive constants $B_{\theta_{e}}$, $B_{q}$,
$C_{\theta_{e}}$, and $C_{q}$ the identities
$\partial_{3}p=B_{\theta_{e}}\theta_{e}-B_{q}qH_{u}$ and
$M=q+\frac{C_{q}}{C_{\theta_{e}}}\theta_{e}$ hold. Then we have the identity
$q=C_{q}\left(\frac{1}{N_{u}^{2}}H_{u}+\frac{1}{N_{s}^{2}}H_{s}\right)\left(\frac{B_{\theta_{e}}C_{\theta_{e}}}{C_{q}}M-\partial_{3}p\right)$
as well as the identities given by

$\displaystyle\frac{\theta_{e}}{C_{\theta_{e}}}$
$\displaystyle=\frac{1}{N_{u}^{2}}\left(\partial_{3}p+B_{q}M\right)H_{u}+\frac{1}{N_{s}^{2}}\left(\partial_{3}p\right)H_{s}$
$\displaystyle=\frac{1}{N_{s}^{2}}\partial_{3}p+\left(\frac{1}{N_{s}^{2}}-\frac{1}{N_{u}^{2}}\right){\min}_{0}\left(\frac{B_{\theta_{e}}C_{\theta_{e}}}{C_{q}}M-\partial_{3}p\right)$

where $N_{u}^{s}$ and $N_{s}^{2}$ are defined in Table 2 and where, as usual,
$\min_{0}(x)\vcentcolon=\min(x,\,0)$.

In particular, since elsewhere in this paper the constants $f$,
$B_{\theta_{e}}$, $B_{q}$, $C_{\theta_{e}}$, and $C_{q}$ are set to unity, we
obtain the following expression for the potential temperature $\theta_{e}$ and
the water content $q$. This expression for $q$ is essential in taking
$PV_{e}$-and-$M$ inversion from its original form (1) to its reformulation (2)
which is used throughout this paper.

######  Corollary 10.

Suppose that $\partial_{3}p=\theta_{e}-qH_{u}$ and $M=\theta_{e}+q$. Then we
have the identities
$\theta_{e}=\partial_{3}p+\frac{1}{2}\min\left(M-\partial_{3}p,\,0\right)$ and
$q=\left(\frac{1}{2}H_{u}+H_{S}\right)\left(M-\partial_{3}p\right)$. In
particular the sign of $q$ agrees with the sign of $M-\partial_{3}p$.

######  Proof C.2.

This follows from setting $B_{\theta_{e}}=B_{q}=C_{\theta_{e}}=C_{q}=1$ in
Lemma C.1.

We conclude this section by noting that the form of $PV_{e}$-and-$M$ inversion
recorded in (1) and studied throughout this paper has nothing to envy to the
more general version (35) recorded above. The latter is _quantitatively_ more
general, but the only _qualitative_ relation that matters between the various
constants at play in (35) is that $N_{u}^{2}>N_{s}^{2}$ (which is always the
case since $N_{u}^{2}-N_{s}^{2}=B_{q}C_{q}>0$). By setting all of the
constants $f$, $B_{\theta_{e}}$, $B_{q}$, $C_{\theta_{e}}$, and $C_{q}$ to
unity (as is done elsewhere in this paper) we are therefore studying the case
where $N_{u}^{2}=2$ and = $N_{s}^{2}=1$. As far as the well-posedness theory
and the iterative methods studied in this paper are concerned, this is
computationally simpler while describing the same qualitative behaviour as any
other choice of values for the constants $f$, $B_{\theta_{e}}$, $B_{q}$,
$C_{\theta_{e}}$, and $C_{q}$.

##  Appendix D Comparison of the variational energy and the conserved energy

In this section we discuss how the variational energy introduced in this paper
in order to characterise $PV_{e}$-and-$M$ inversion (c.f. Definition 4.1)
differs from the energy that is conserved by the dynamics of the precipitating
quasi-geostrophic equations.

Indeed, while this conserved energy may at first glance seem like a viable
candidate for the variational formulation of $PV_{e}$-and-$M$ inversion, we
will see that it is not the case. We will also see that the energy we
introduce in Definition 4.1 has better regularity properties than the
conserved energy. First we write the conserved energy in terms of solely $p$
and $M$.

######  Lemma D.1 (Rewriting the conserved energy).

The conserved energy $E_{cons}$ introduced in [MSS19] as $E_{1}$ may, in the
notation of Appendix C, be written as follows:

$E_{cons}=\int_{\mathbb{T}^{3}}\frac{1}{2}{\lvert
u\rvert}^{2}+\left(\frac{\widetilde{B}_{\theta_{e}}}{2N_{u}^{2}}H_{u}+\frac{B_{\theta_{e}}}{2N_{s}^{2}}H_{s}\right)\theta_{e}^{2}+\frac{\widetilde{B}_{q}}{2N_{s}^{2}}q^{2}H_{u}$
(36)

for
$\widetilde{B}_{\theta_{e}}=B_{\theta_{e}}\left(B_{\theta_{e}}+B_{q}C_{q}/C_{\theta_{e}}\right)$
and
$\widetilde{B}_{q}=B_{q}\left(B_{q}+B_{\theta_{e}}C_{\theta_{e}}/C_{q}\right)$.
In particular, if we consider $B_{\theta_{e}}=B_{q}=C_{\theta_{e}}=C_{q}=1$
then $E_{cons}=\int_{\mathbb{T}^{3}}\frac{1}{2}{\lvert
u\rvert}^{2}+\frac{1}{2}\theta_{e}^{2}+\frac{1}{2}q^{2}H_{u}$ such that, if
geostrophic balance $u=\nabla_{h}^{\perp}p$ and hydrostatic balance
$\partial_{3}p={\theta_{e}}-qH_{u}$ hold,

$E_{cons}=\int_{\mathbb{T}^{3}}\frac{1}{2}{\lvert\nabla
p\rvert}^{2}+\frac{1}{4}(M+\partial_{3}p){\min}_{0}\,(M-\partial_{3}p)$ (37)

where as usual, $\min_{0}(x)\vcentcolon=\min(x,\,0)$.

######  Proof D.2.

We recall from [MSS19] that the conserved energy is defined to be

$E_{cons}\vcentcolon=\int_{\mathbb{T}^{3}}\frac{1}{2}{\lvert
u\rvert}^{2}+\frac{b_{u}^{2}}{2N_{u}^{2}}H_{u}+\frac{b_{s}^{2}}{2N_{s}^{2}}H_{s}+\frac{N_{u}^{2}N_{s}^{2}}{2(N_{u}^{2}-N_{s}^{2})}\widetilde{M}^{2}H_{u}$

where $b_{u}=\theta_{e}-q$ and $b_{s}=\theta_{e}$ such that the buoyancy may
be decomposed as $\theta_{e}-qH_{u}=b_{u}H_{u}+b_{s}H_{s}$ and where
$\widetilde{M}$ is a constant multiple of the moist variable $M$ used herein.
Indeed we have that
$\widetilde{M}\vcentcolon=b_{u}/N_{u}^{2}-b_{s}/N_{s}^{2}=-(B_{q}/N_{u}^{2})M.$
A straightforward but tedious computation then produces (36). To obtain (37)
we use Lemma C.1 to write $\theta_{e}$ and $q$ in terms of $M$ and
$\partial_{3}p$ and plug into the expression for the conserved energy. A short
computation then verifies that

$\displaystyle{\lvert
u\rvert}^{2}+\theta_{e}^{2}+q^{2}H_{u}={\lvert\nabla_{h}p\rvert}^{2}+{\left(\partial_{3}p+\frac{1}{2}(M-\partial_{3}p)H_{u}\right)}^{2}+\frac{1}{4}{(M-\partial_{3}p)}^{2}H_{u}$
$\displaystyle={\lvert\nabla
p\rvert}^{2}+\frac{1}{2}(M+\partial_{3}p){\min}_{0}\,(M-\partial_{3}p),$

as desired.

We may then use (37) to compute the Euler-Lagrange equation associated with
the conserved energy.

######  Corollary 11 (Euler-Lagrange equation associated with the conserved
energy).

Let $E_{cons}$ be as defined as in Lemma D.1, with $M$ being fixed and viewing
$E_{cons}$ as a function of the pressure $p$ only. The Euler-Lagrange equation
associated with $E_{cons}$ is $-\Delta
p+\frac{1}{2}\partial_{3}((\partial_{3}p)H_{u})=0.$

######  Proof D.3.

We simply take a derivative of $E_{cons}(p)$ in the direction of $\phi$,
obtaining $DE_{cons}(p)\phi=\int_{\mathbb{T}^{3}}\nabla
p\cdot\nabla\phi-\frac{1}{2}(\partial_{3}p)H_{u}(\partial_{3}\phi),$ from
which the claim follows.

######  Remark D.4 (Comparison of the variational energy and the conserved
energy).

Combining Lemma 4.7 and Corollary 11 tells us that

$DE(p)\phi=\underbrace{DE_{cons}(p)\phi}_{\sim\,\text{principal
part}}+\underbrace{\langle
PV_{e}-\frac{1}{2}\partial_{3}(MH_{u}),\,p\rangle}_{\text{nonlinear
forcing}}.$ (38)

What is very interesting about this identity is that, in some sense, the
conserved energy extracts the principal part of the operator characterising
$PV_{e}$-and-$M$ inversion. Indeed, we can essentially view
$PV_{e}-\frac{1}{2}\partial_{3}(MH_{u})$ as a term forcing the Euler-Lagrange
operator associated with the conserved energy. However, and this is crucial,
this forcing term is _nonlinear_ in $p$ due to the presence of $H_{u}$. The
derivative of $\frac{1}{2}\langle\partial_{3}(MH_{u}),\,p\rangle$ with respect
to $p$ is _not_ equal to $\frac{1}{2}\partial_{3}(MH_{u})$ itself which means
that we cannot simply claim that $PV_{e}$-and-$M$ inversion is the Euler-
Lagrange equation associated with
$E_{cons}+PV_{e}-\frac{1}{2}\partial_{3}(MH_{u})$.

This means that, in order to identify a variational formulation of
$PV_{e}$-and-$M$ inversion it is not sufficient to take inspiration from the
conserved energy. We must instead analyse the PDE carefully, which leads to
the definition of a _different_ energy. The identity (38) also highlights the
essential trick underlying the definition of the variational energy $E$: the
forcing term $\frac{1}{2}\partial_{3}(MH_{u})$ is carefully folded into the
energy, thus obtaining a valid variational formulation of $PV_{e}$-and-$M$
inversion.

In particular this leads to a more regular energy. Indeed a comparison of (4)
and (37) shows that the energy density appearing in the variational energy is
$C^{1,\,1}$, as a function of $\partial_{3}p$, whereas the energy density
appearing in the conserved energy is only $C^{0,\,1}$.

## Acknowledgments

This work was partially supported by the National Science Foundation under
grant NSF-DMS-1907667. As recipient of an INI-Simons Postdoctoral Research
Fellowship, A.R.-T. gratefully acknowledges support from the Simons
Foundation, the Isaac Newton Institute for Mathematical Sciences, and the
Department of Applied Mathematics and Theoretical Physics of the University of
Cambridge.

The authors thank Peter Constantin for helpful discussions in the early stages
of this work.

## References

  * [BCZT14] A. Bousquet, M. Coti Zelati, and R. Temam.  Phase transition models in atmospheric dynamics.  Milan J. Math., 82(1):99–128, 2014.
  * [Cha48] J. G. Charney.  On the scale of atmospheric motions.  Geofys. Publ. Norske Vid.-Akad. Oslo, 17(2):17, 1948.
  * [CHT+18] Y. Cao, M. Hamouda, R. Temam, J. Tribbia, and X. Wang.  The equations of the multi-phase humid atmosphere expressed as a quasi variational inequality.  Nonlinearity, 31(10):4692–4723, 2018.
  * [CJTT21] Y. Cao, C. Jia, R. Temam, and J. Tribbia.  Mathematical analysis of a cloud resolving model including the ice microphysics.  Discrete Contin. Dyn. Syst., 41(1):131–167, 2021.
  * [CSV15] G.-Q. Chen, H. Shahgholian, and J.-L. Vazquez.  Free boundary problems: the forefront of current and future developments.  Trans. R. Soc. A, 373, 2015.
  * [CT07] C. Cao and E. S. Titi.  Global well-posedness of the three-dimensional viscous primitive equations of large scale ocean and atmosphere dynamics.  Ann. of Math. (2), 166(1):245–267, 2007.
  * [CV10] L. A. Caffarelli and A. F. Vasseur.  The De Giorgi method for regularity of solutions of elliptic equations and its applications to fluid dynamics.  Discrete Contin. Dyn. Syst. Ser. S, 3(3):409–427, 2010.
  * [CZFTT13] M. Coti Zelati, M. Frémond, R. Temam, and J. Tribbia.  The equations of the atmosphere with humidity and saturation: uniqueness and physical bounds.  Phys. D, 264:49–65, 2013.
  * [CZHK+15] M. Coti Zelati, A. Huang, I. Kukavica, R. Temam, and M. Ziane.  The primitive equations of the atmosphere in presence of vapour saturation.  Nonlinearity, 28(3):625–668, 2015.
  * [CZT12] M. Coti Zelati and R. Temam.  The atmospheric equation of water vapor with saturation.  Boll. Unione Mat. Ital. (9), 5(2):309–336, 2012.
  * [DG57] E. De Giorgi.  Sulla differenziabilita e l’analiticita della estremali degli integrali multipli regolari.  Mem. Accad. Sci. Torino. Cl. Sci. Fis. Mat. Nat., Vol. 3:25–43, 1957.
  * [Eck60] C. Eckart.  Hydrodynamics of oceans and atmospheres.  Pergamon Press, New York-Oxford-London-Paris, 1960.
  * [Ert42] H. Ertel.  Ein neuer hydrodynamischer wirbelsatz.  Meteor. Zeitschrift, 59:277–281, 1942.
  * [ESS20a] T. K. Edwards, L. M. Smith, and S. N. Stechmann.  Atmospheric rivers and water fluxes in precipitating quasi-geostrophic turbulence.  Quarterly Journal of the Royal Meteorological Society, 146(729):1960–1975, 2020.
  * [ESS20b] T. K. Edwards, L. M. Smith, and S. N. Stechmann.  Spectra of atmospheric water in precipitating quasi-geostrophic turbulence.  Geophys. Astrophys. Fluid Dyn., 114(6):715–741, 2020.
  * [Fol99] G. B. Folland.  Real analysis.  Pure and Applied Mathematics (New York). John Wiley & Sons, Inc., New York, second edition, 1999.  Modern techniques and their applications, A Wiley-Interscience Publication.
  * [FS15] A. Figalli and H. Shahgholian.  An overview of unconstrained free boundary problems.  Trans. R. Soc. A, 373, 2015.
  * [GM12] M. Giaquinta and L. Martinazzi.  An introduction to the regularity theory for elliptic systems, harmonic maps and minimal graphs, volume 11 of Appunti. Scuola Normale Superiore di Pisa (Nuova Serie) [Lecture Notes. Scuola Normale Superiore di Pisa (New Series)].  Edizioni della Normale, Pisa, second edition, 2012.
  * [Gra98] W. W. Grabowski.  Toward cloud resolving modeling of large-scale tropical circulations: A simple cloud microphysics parametrisation.  Journal of Atmospheric Sciences, 55:3283–3298, 1998.
  * [GT01] D. Gilbarg and N. S. Trudinger.  Elliptic partial differential equations of second order.  Classics in Mathematics. Springer-Verlag, Berlin, 2001.  Reprint of the 1998 edition.
  * [HDMSS13] G. Hernandez-Duenas, A. J. Majda, L. M. Smith, and S. N. Stechmann.  Minimal models for precipitating turbulent convection.  J. Fluid Mech., 717:576–611, 2013.
  * [HESS21] R. Hu, T. K. Edwards, L. M. Smith, and S. N. Stechmann.  Initial investigations of precipitating quasi-geostrophic turbulence with phase changes.  Res. Math. Sci., 8(1):Paper No. 6, 25, 2021.
  * [HKLT17] S. Hittmeir, R. Klein, J. Li, and E. S. Titi.  Global well-posedness for passively transported nonlinear moisture dynamics with phase changes.  Nonlinearity, 30(10):3676–3718, 2017.
  * [HKLT20] S. Hittmeir, R. Klein, J. Li, and E. S. Titi.  Global well-posedness for the primitive equations coupled to nonlinear moisture dynamics with phase changes.  Nonlinearity, 33(7):3206–3236, 2020.
  * [HL11] Q. Han and F. Lin.  Elliptic partial differential equations, volume 1 of Courant Lecture Notes in Mathematics.  Courant Institute of Mathematical Sciences, New York; American Mathematical Society, Providence, RI, second edition, 2011.
  * [HMR85] B. J. Hoskins, M. E. McIntyre, and A. W. Robertson.  On the use and significance of isentropic potential vorticity maps.  Quarterly Journal of the Royal Meteorological Society, 111(470):877–946, October 1985.
  * [Leo17] G. Leoni.  A first course in Sobolev spaces, volume 181 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, second edition, 2017.
  * [LM20] R. Lian and J. Ma.  Existence of a strong solution to moist atmospheric equations with the effects of topography.  Bound. Value Probl., pages Paper No. 103, 34, 2020.
  * [Maj03] A. Majda.  Introduction to PDEs and waves for the atmosphere and ocean, volume 9 of Courant Lecture Notes in Mathematics.  New York University, Courant Institute of Mathematical Sciences, New York; American Mathematical Society, Providence, RI, 2003.
  * [MK06] A. J. Majda and R. Klein.  Systematic multiscale models for deep convection on mesoscales.  Theor. Comput. Fluid Dyn., 20:525–551, 2006.
  * [MSS19] D. H. Marsico, L. M. Smith, and S. N. Stechmann.  Energy decompositions for moist boussinesq and anelastic equations with phase changes.  Journal of the Atmospheric Sciences, 76(11):3569 – 3587, 2019.
  * [Nes04] Y. Nesterov.  Introductory lectures on convex optimization, volume 87 of Applied Optimization.  Kluwer Academic Publishers, Boston, MA, 2004.  A basic course.
  * [PSU12] A. Petrosyan, H. Shahgholian, and N. Uraltseva.  Regularity of free boundaries in obstacle-type problems, volume 136 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, 2012.
  * [SS17] L. M. Smith and S. N. Stechmann.  Precipitating quasigeostrophic equations and potential vorticity inversion with phase changes.  Journal of the Atmospheric Sciences, 74(10):3285 – 3303, 2017.
  * [TL22] S. Tan and W. Liu.  The strong solutions to the primitive equations coupled with multi-phase moisture atmosphere.  Phys. D, 440:Paper No. 133442, 18, 2022.
  * [TW15] R. Temam and K. J. Wu.  Formulation of the equations of the humid atmosphere in the context of variational inequalities.  J. Funct. Anal., 269(7):2187–2221, 2015.