Version of August 28, 2024$,\,\,\,$

# A preliminary model for optimal control of moisture content in unsaturated
soils

Marco Berardi  Istituto di Ricerca sulle Acque, Consiglio Nazionale delle
Ricerche, Via De Blasio, 5, 70132 Bari, Italy<EMAIL_ADDRESS>,
Fabio V. Difonzo  Dipartimento di Matematica, Università degli Studi di Bari
Aldo Moro, Via E. Orabona 4, 70125 Bari, Italy<EMAIL_ADDRESS>and
Roberto Guglielmi  Department of Applied Mathematics, University of Waterloo,
200 University Ave W, Waterloo, N2L 3G1 Ontario, Canada
<EMAIL_ADDRESS>
###### Abstract.

In this paper we introduce an optimal control approach to Richards’ equation
in an irrigation framework, aimed at minimizing water consumption while
maximizing root water uptake. We first describe the physics of the nonlinear
model under consideration, and then develop the first-order necessary
optimality conditions of the associated boundary control problem. We show that
our model provides a promising framework to support optimized irrigation
strategies, thus facing water scarcity in irrigation. The characterization of
the optimal control in terms of a suitable relation with the adjoint state of
the optimality conditions is then used to develop numerical simulations on
different hydrological settings, that supports the analytical findings of the
paper.

###### Key words and phrases:

Richards’ equation, optimal control in agriculture, direct-dual optimality
system, optimal control of boundary conditions, modeling and control in soil
systems

###### 1991 Mathematics Subject Classification:

34H05, 76S05

##  1\. Introduction

More and more often, extreme weather events are accompanied by longer, more
intense heat waves and consequent periods of drought, and a forecast global
warming is increasing the urgent need of freshwater for human life. In this
context, freshwater necessary for agriculture represents almost 70% of the
whole amount of freshwater reserve [1].

In this scenario, a wise management of water resources for agricultural
purposes is of fundamental importance, even at the irrigation district scale
[2]. Nevertheless, the vast majority of irrigation models just applies
heuristic approaches for determining the amount and the timing of irrigation.
Albeit an expert knowledge of agricultural and phenological issue is crucial,
very seldom such information is coupled with proper mathematical models for
controlling irrigation. In most sophisticated cases, the irrigation is managed
mainly by studying soil water infiltration into the root zone, starting from
Gardner’s pioneering works, reviewed in [3]; several tools have been proposed
in this context, generally providing some type of solver for Richards’
equation, the advection-diffusion equation which describes water infiltration
in unsaturated porous media accounting also for root water uptake models.
In the water resources management or agronomic framework several tools have
been proposed in order to benefit from the Richards’ equation for irrigation
purposes: in [4] a Python code is presented, able to solve the Richards’
equation with any type of root uptake model by a transverse method of lines;
also, the popular Hydrus software is often used for simulating water flow and
root uptake with different crops (e.g apples, as in [5], pecan trees as in
[6]).

To the best of our knowledge, control methods are seldom applied to irrigation
problems, and always with simplified models. For instance, in [1] a zone model
predictive control is designed after defining a linear parameter varying
model, aimed at maintaining the soil moisture in the root zone within a
certain target interval; in [7] an optimal control is applied to an irrigation
problem, modelled by a simpler (with respect to Richards’ equation) hydrologic
balance law; a simplified optimization method based on the computation of
steady solutions of Richards’ equation is proposed in [8]; finally, [9]
presents an interesting sliding-mode control approach, but only considering a
constant diffusion term in Richards’ equation. An elegant approach for
applying control techniques in a Richards’ equation framework is provided in
[10], yet with very different applications and tools, i.e. maximizing the
amount of absorbed liquid by redistributing the materials, when designing the
material properties of a diaper.

In this paper, we propose a model for solving an optimal control problem under
the quasi-unsaturated assumptions (see Section 2), which provide a suitable
hydrological setting that prevents to reach water moisture saturation in the
soil. We derive the appropriate optimality conditions for the boundary control
of a class of nonlinear Richards’ equations, and implement these results in
the development and computation of numerical solutions by a classical
Projected Gradient Descent algorithm.
Albeit the focus of this paper does not consist in proposing a novel numerical
method, and here a standard MATLAB solver is integrated with control tools,
some significant advances in the numerical solution of the unsaturated flow
model deserve to be reminded; as a matter of fact, the numerics literature on
Richards’ equation is currently enriching and constantly evolving, since its
possibly degenerate and highly nonlinear nature poses several challenges. For
instance, the treatment of nonlinearities is a significant issue, and has been
faced by different techniques, as Newton methods ([11, 12]), L-scheme or its
variants [13, 14, 15] or Picard iterations [16]. The discretization in space
has been dealt, for instance, by finite elements or mixed finite element
methods [17, 18, 19], discontinuous Galerkin [20, 21], finite volume methods
[22, 23]. A separate mention is deserved by the problem of infiltration in
presence of discontinuities, which can be handled by domain decomposition
methods [24, 25], Filippov approach [26, 27]. For a more detailed discussion,
the interested reader is addressed to the following complete reviews on
numerical issues in Richards’ equation [28, 29]. The aforementioned methods
can be used to face peculiar issues in the numerical integration of Richards’
equation, where MATLAB pdepe is known to blow up in a relatively small
integration time (see, e.g., [4]). Thus, they could provide further directions
in the development of specific algorithms for solving optimal control
problems.

It is worth stressing that there is a vast literature about the problems of
existence, uniqueness and regularity of the solution to degenerate parabolic
differential equations. In the specific case of Richards’ equations, the
existence of solutions is tackled in the seminal paper by Alt and Luckhaus
[30], whose ideas were subsequently developed by other authors. For example,
[31] describes the semigroup approach to determine the existence of weak
solutions, further developed in [32]. In this paper we adopt the functional
framework from [33], that exhaustively describes the minimal assumptions to
derive the necessary regularity conditions to develop our analysis, for
general classes of hydrological settings. We also refer to [34] for a thorough
review of such results.

The paper is organized as follows: In Section 2 we introduce the quasi-
unsaturated model of the Richards’ equation; in Section 3 we describe the
framework to ensure the well-posedness of the model of interest, and then
derive first order necessary optimality conditions via the Lagrangian method
in Section 4. Finally, we present numerical simulations in Section 5.

##  2\. The mathematical model

Our work stems from the quasi-unsaturated Richards’ model [33], describing a
fast diffusion of water in soils. This framework provides a convenient setting
to apply optimal control methods and derive optimal irrigation strategies.
Indeed, from the mathematical point of view, the quasi-unsaturated diffusive
model retains many crucial features of the nonlinear diffusion of water in the
soil, while avoiding the special mathematical treatment required to face the
limit case of a saturated diffusion. We consider the model of water diffusion
in the space domain $(0,Z)$, where $Z>0$ is the depth of the domain under
consideration. We denote by $T\in(0,+\infty)$ the time horizon of the interval
$(0,T)$, by $Q=(0,Z)\times(0,T)$ the space-time domain, and by $(\cdot,\cdot)$
and $\|\,\cdot\,\|$ the scalar product and the norm in $L^{2}(\Omega)$,
respectively. In terms of hydraulic parameters,
$\theta:Q\to[\theta_{r},\theta_{S})$ is the _water content_ or _moisture_ ,
where $\theta_{r}$ and $\theta_{S}$ represent the _residual_ and the
_saturated water content_ , respectively. The function
$\beta:[\theta_{r},\theta_{S})\to[\varrho,+\infty)$ is the _water diffusivity_
, satisfying the following condition:

  * $\mathbf{(H_{\beta})}$

$\beta$ is locally Lipschitz continuous and monotonically increasing,
$\beta(\theta)\geq\varrho>0$ for all $\theta\in[\theta_{r},\theta_{S})$, and
$\displaystyle\lim_{\theta\nearrow\theta_{S}}\beta(\theta)=+\infty$.

The function $\beta^{*}$ is the primitive of the water diffusivity $\beta$
that vanishes at $0$. Thus, assumption $(\mathbf{H}_{\beta})$ implies that
$\beta^{*}$ is differentiable and monotonically increasing on
$[\theta_{r},\theta_{S})$, and satisfies

$(\beta^{*}(\theta_{1})-\beta^{*}(\theta_{2}))(\theta_{1}-\theta_{2})\geq\varrho(\theta_{1}-\theta_{2})^{2}\qquad\forall\;\theta_{1},\theta_{2}\in[\theta_{r},\theta_{S})\,.$

Moreover, the _hydraulic conductivity_
$K:[\theta_{r},\theta_{S}]\to\mathbb{R}$ is a non-negative, Lipschitz
continuous on $[\theta_{r},\theta_{S}]$, and monotonically increasing
function. Other hydraulic functions of interest are the _liquid pressure head
$h:Q\to(-\infty,0)$_, that is negative for unsaturated porous media, and the
_specific water capacity_ $C(h)=\frac{\mathrm{d}\theta}{\mathrm{d}h}$, that
practically represents a storage term. The relation between the functions
$\beta$, $K$ and $C$ is then expressed by
$\beta(\theta(h)):=\frac{K(h)}{C(h)}$.
With these notations, and assuming the vertical axis with downward positive
orientation, the implicit form of the quasi-unsaturated model of the
nonhysteretic infiltration of an incompressible fluid into an isotropic,
homogeneous, unsaturated porous medium with a constant porosity and truncated
diffusivity with non-homogeneous–Dirichlet Boundary Conditions (BCs) is given
by the system

(2.1) $\begin{cases}\frac{\partial\theta}{\partial
t}-\frac{\partial^{2}\beta^{*}(\theta)}{\partial z^{2}}+\frac{\partial
K(\theta)}{\partial z}=f(\theta)&\text{ in }\ Q\,,\\\
\theta(z,0)=\theta_{0}(z)&\text{ in }\ (0,Z)\,,\\\ \theta(0,t)=v(t)&\text{ for
}\ t\in(0,T)\,,\\\ \theta(Z,t)=g(t)&\text{ for }\ t\in(0,T)\,.\end{cases}$

The source term $f(\theta)$ represents a sink function, that in our model
describes the root water uptake. The non-homogeneous Dirichlet BCs are given
by two functions $v,g:(0,T)\to\mathbb{R}$, where the BCs at $z=0$ is the
control input $v$, that describes the irrigation strategy over the time
horizon $(0,T)$, while the BCs at $z=Z$ is a given function $g$ modeling the
interaction with the environment below the root zone

(2.2) $\theta_{r}<v(t),g(t)<\theta_{S}\quad\text{ for a.e. }\ t\in(0,T)\;.$

######  Remark 2.1.

Let us notice that the first equation in (2.1) is the Richards’ equation. In
fact, we can easily compute that

$\displaystyle\frac{\partial}{\partial z}\beta^{*}$
$\displaystyle=\frac{\partial\beta^{*}}{\partial\theta}\frac{\partial\theta}{\partial
h}\frac{\partial h}{\partial z}=\frac{K(h)}{C(h)}C(h)\frac{\partial
h}{\partial z}=K(h)\frac{\partial h}{\partial z},$

that provides the diffusion term in the classical mixed form of the Richards’
equation [16, Eq. (3)].

##  3\. Existence of solutions

Following [33, Section 4.5], we first reduce system (2.1) to a problem with
homogeneous Dirichlet BCs. For this purpose, we assume the following
hypothesis

$\displaystyle\mathbf{(H_{\omega})}\quad\exists\;\omega:Q\to\mathbb{R}\ \text{
such that }$ $\displaystyle\begin{cases}\omega\in L^{2}(0,T;H^{1}(0,Z))\cap
L^{\infty}(Q)\,,\quad\omega_{t}\in L^{2}(Q)\,,&\\\
\|\omega\|_{L^{\infty}(Q)}<\theta_{S}\,,\quad\mathrm{essinf}_{(z,t)\in
Q}\;\omega(z,t)>\theta_{r}\;,&\\\ \omega(0,t)=v(t)\,,\ \omega(Z,t)=g(t)\,,\
\text{ for a.e. }\ t\in(0,T)\,,&\end{cases}$

where $\omega_{t}$ denotes the time derivative in the sense of distributions
from $(0,T)$ to $L^{2}(0,Z)$. Thus, the function $\omega$ is defined on the
cylinder $Q$ and it attains the boundary values of $v$ and $g$ in $z=0$ and
$z=Z$, respectively. Then, introducing the function $\phi=\theta-\omega$,
system (2.1) is equivalent to

(3.1) $\left\\{{\begin{array}[]{ll}\frac{\partial\phi}{\partial
t}-\frac{\partial^{2}F^{\omega}(\phi)}{\partial z^{2}}+\frac{\partial
K(\phi+\omega)}{\partial z}=f^{B}(\phi)-\omega_{t}&\text{ in }\
Q\,,\\\\[3.44444pt] \phi(z,0)=\phi_{0}(z)&\text{ in }\ (0,Z)\,,\\\\[1.72218pt]
\phi(0,t)=0&\ t\in(0,T)\,,\\\\[1.72218pt] \phi(Z,t)=0&\
t\in(0,T)\,,\end{array}}\right.$

where $\phi_{0}:=\theta_{0}-\omega_{0}$ and, for all $\phi\in V$,

$F^{\omega}(\phi):=\beta^{*}(\phi+\omega)-\beta^{*}(\omega)\;,\qquad
f^{B}(\phi):=f(\phi+\omega)+\frac{\partial^{2}\beta^{*}(\omega)}{\partial
z^{2}}\;.$

We now introduce a suitable functional framework for (3.1). Let us consider
the Gelfand triple $H=L^{2}(0,Z)$, $V=H^{1}_{0}(0,Z)$, and its dual
$V^{\prime}=H^{-1}(0,Z)$, with their usual norms. Then, (3.1) is equivalent to
the abstract equation

(3.2)
$\frac{\mathrm{d}\phi}{\mathrm{d}t}+B(t)\phi=f^{B}(\phi)-\omega_{t}\quad\text{
a.e. }\ t\in(0,T)\,,\qquad\phi(0)=\phi_{0}\,,$

where the operator $B(t):V\to V^{\prime}$ is defined by

$\langle B(t)\phi,\psi\rangle_{V^{\prime},V}=\int_{0}^{Z}\left(\frac{\partial
F^{\omega}(\phi(t))}{\partial
z}-K(\phi(t)+\omega(t))\right)\frac{\partial\psi}{\partial
z}\mathrm{d}z\qquad\forall\phi,\psi\in V\;.$

Notice that, thanks to $\mathbf{(H_{\omega})}$, we have that
$\frac{\partial^{2}\beta^{*}(\omega)}{\partial z^{2}}\in
L^{2}(0,T;V^{\prime})$. Indeed, for all $\varphi\in V=H^{1}_{0}(0,Z)$,

$\left\lvert\int_{0}^{Z}-\frac{\partial^{2}\beta^{*}(\omega)}{\partial
z^{2}}\varphi\,\mathrm{d}z\right\rvert=\left\lvert\int_{0}^{Z}\frac{\partial\beta^{*}(\omega)}{\partial
z}\frac{\partial\varphi}{\partial
z}\mathrm{d}z\right\rvert\leq\left\|\beta(\omega)\frac{\partial\omega}{\partial
z}\right\|\,\left\|\frac{\partial\varphi}{\partial z}\right\|\leq
M_{\omega}\|\varphi\|_{V}\,.$

Hence $\|\frac{\partial^{2}\beta^{*}(\omega)}{\partial
z^{2}}\|_{V^{\prime}}\leq M_{\omega}$, and thus $f^{B}-\omega_{t}\in
L^{2}(0,T;V^{\prime})$. In particular, this implies that the right-hand side
of (3.2) is also in $L^{2}(0,T;V^{\prime})$.

######  Definition 3.1.

Let $\theta_{0}\in L^{2}(0,Z)$ and $f:[\theta_{r},\theta_{S}]\to\mathbb{R}$ be
Lipschitz continuous. We say that the function $\phi\in C([0,T];L^{2}(0,Z))$
is a solution to (3.2) if $\frac{\mathrm{d}\phi}{\mathrm{d}t}\in
L^{2}(0,T;V^{\prime})$, $F^{\omega}\in L^{2}(0,T;V)$, $\phi(0)=\phi_{0}$ and,
for a.e. $t\in(0,T)$ and for all $\psi\in V$,

(3.3)
$\left\langle\frac{\mathrm{d}\phi}{\mathrm{d}t}(t),\psi\right\rangle_{V^{\prime},V}+\langle
B(t)\phi,\psi\rangle_{V^{\prime},V}=\left\langle
f^{B}(\phi)-\omega_{t},\psi\right\rangle_{V^{\prime},V}\;.$

Existence of solutions to (3.2) with a source term independent of the water
content $\theta$ is proved in [33, Proposition 5.3]. For our purpose, we need
to extend such well-posedness result to the case of a nonlinear source term
$f(\theta)$, as in system (2.1). This can indeed be achieved by following a
standard Galerkin approximation approach (see, e.g., [35, Lemma 5.3]).

######  Theorem 3.2.

Assume $\mathbf{(H_{\omega})}$, $\theta_{0}\in L^{2}(0,Z)$, and that
$f:[\theta_{r},\theta_{S}]\to\mathbb{R}$ is Lipschitz continuous. Then the
problem (3.2) admits a unique solution $\phi\in C([0,T];L^{2}(0,Z))\cap
L^{2}(0,T;V)$ with $\frac{\mathrm{d}\phi}{\mathrm{d}t}\in
L^{2}(0,T;V^{\prime})$ and $F^{\omega}(\phi)\in L^{2}(0,T;V)$. Therefore,
system (2.1) admits a unique solution $\theta\in L^{2}(0,T;H^{1}(0,Z))\cap
C([0,T];H)$ with $\frac{\mathrm{d}\theta}{\mathrm{d}t}\in
L^{2}(0,T;V^{\prime})$ and $\beta^{*}(\theta)\in L^{2}(0,T;H^{1}(0,Z))$.

Moreover, for appropriate initial conditions, we can prove that the solution
stays away from the saturation value $\theta_{S}$ uniformly in time. To this
aim, we introduce the function $j:\mathbb{R}\to(-\infty,\infty]$ defined by

$j(r)\mathrel{\mathop{:}}=\begin{cases}\displaystyle\int_{0}^{r}\beta^{*}(\xi)d\xi\,,&r<\theta_{S}\,,\\\\[2.15277pt]
\displaystyle\infty\,,&r\geq\theta_{S}\,,\end{cases}$

and the space

$M_{j}\mathrel{\mathop{:}}=\\{\theta\in L^{2}(0,Z):j(\theta)\in
L^{1}(0,Z)\\}\,.$

######  Theorem 3.3 (Theorem 5.7, [33]).

Assume $\mathbf{(H_{\omega})}$, $\theta_{0}\in M_{j}$, and that
$f:[\theta_{r},\theta_{S}]\to\mathbb{R}$ is Lipschitz continuous. Moreover,
assume that $f$ is non-negative, that is, there exists $f_{m}\in[0,\infty)$
such that $f_{m}\leq f$, and

$\displaystyle\mathrm{essinf}_{x\in(0,Z)}\theta_{0}(x)\geq 0\,,$
$\displaystyle\theta_{m}(t)\leq g(t),v(t)<\theta_{S}\,,\quad\text{for all }\
t\in[0,T]\,,$

where

$\theta_{m}(t)=\mathrm{essinf}_{x\in(0,Z)}\theta_{0}(x)+f_{m}t\,.$

Then the solution $\theta$ to problem (2.1) satisfies

$\theta_{m}(t)\leq\theta(x,t)<\theta_{S}\,,\quad\text{for all }\
(x,t)\in\\{0,Z\\}\times[0,T]\,.$

##  4\. The optimal control problem

In this section we formally derive the first order necessary optimality
conditions for the cost functional

(4.1) $J(\theta,u)=\frac{1}{2}\int_{Q}\left\lvert
f(\theta(z,t))-1\right\rvert^{2}\,\mathrm{d}z\mathrm{d}t+\frac{\lambda}{2}\int_{0}^{T}\left\lvert
u(t)\right\rvert^{2}\,\mathrm{d}t\;,$

where $u=v-\theta_{r}$, $v$ is the control that appears in (2.1), $\lambda>0$
is the coefficient of the control cost,
$f:[\theta_{r},\theta_{S}]\to\mathbb{R}$ describes the normalized root water
uptake model as in (2.1), and $\theta$ is the solution to (2.1) with $f$ as
the sink term. Roughly speaking, the performance index (4.1) optimizes the
root water uptake (see, for example, expression (5.2) in Section 5, where $f$
is maximized when $f\equiv 1$) while minimizing the irrigation cost $u$. In
this setting, it is natural to consider the following space of admissible
control

(4.2) $U_{ad}:=\\{u\in L^{\infty}(0,T):0\leq u(t)<\theta_{S}-\theta_{r}\
\text{ for a.e. }\ t\in(0,T)\\}\;.$

Fixing $g\in L^{2}(0,T)$, $\theta_{0}\in L^{2}(0,Z)$, we introduce the
control-to-state operator $\Lambda:U_{ad}\to C([0,T];L^{2}(0,Z))$ such that
$u\in U_{ad}\mapsto\theta\in C([0,T];L^{2}(0,Z))$ solution of (2.1). Theorem
3.2 ensures that the mapping $\Lambda$ is well-posed. We can thus reformulate
the minimization of a functional $\widetilde{J}(\theta,u)$ constrained to the
control system (2.1) in terms of the so-called _reduced cost functional_
$J:U_{ad}\to\mathbb{R}$ defined by $J(u):=\widetilde{J}(\Lambda(u),u)$. We
first introduce the Lagrangian functional

$\displaystyle\mathcal{L}(\theta,u,{\bf p})$
$\displaystyle=J(\theta,u)-\int_{Q}\left[\frac{\partial\theta}{\partial
t}-\frac{\partial}{\partial z}\left(\beta\frac{\partial\theta}{\partial
z}\right)+\frac{\partial K(\theta)}{\partial
z}-f\right]p\,\mathrm{d}z\mathrm{d}t$
$\displaystyle\quad-\int_{0}^{T}\left(\theta(0,t)-u(t)\right)p_{1}\,\mathrm{d}t-\int_{0}^{T}\left(\theta(Z,t)-g(t)\right)p_{2}\,\mathrm{d}t,$

where ${\bf p}=(p,p_{1},p_{2})$ are adjoint variables that will be useful to
find a representation of the optimal control. After integration by parts, we
can rewrite the Lagrangian functional as

$\displaystyle\mathcal{L}(\theta,u,{\bf p})$
$\displaystyle=J(\theta,u)-\int_{Q}\frac{\partial\theta}{\partial
t}p+\beta\frac{\partial\theta}{\partial z}\frac{\partial p}{\partial
z}+\left(\frac{\partial K(\theta)}{\partial
z}-f\right)p\,\mathrm{d}z\mathrm{d}t$
$\displaystyle\quad+\int_{0}^{T}\left[\left(\beta\frac{\partial\theta}{\partial
z}p\right)_{\big{|}z=Z}-\left(\beta\frac{\partial\theta}{\partial
z}p\right)_{\big{|}z=0}\right]\,\mathrm{d}t+$
$\displaystyle\quad-\int_{0}^{T}\left(\theta(0,t)-u(t)\right)p_{1}\,\mathrm{d}t-\int_{0}^{T}\left(\theta(Z,t)-g(t)\right)p_{2}\,\mathrm{d}t\;.$

Hereafter, we shall assume that the source term $f\in
H^{1}(\theta_{r},\theta_{S})$ to justify the following computations. In order
to derive the first order optimality conditions of problem (4.1)-(2.1) with
input constraints (4.2), we enforce the condition
$D_{\theta}\mathcal{L}(\theta^{*},u^{*},p^{*})\theta=0$ for all $\theta$, that
determines the equation satisfied by the adjoint variable $p$; and the
condition $D_{u}\mathcal{L}(\theta^{*},u^{*},p^{*})\cdot(u-u^{*})\geq 0$ for
all $u\in U_{ad}$, that returns the optimality condition satisfied by any
optimal control $u^{*}$. After direct computations, we get that

$D_{\theta}\mathcal{L}(\theta^{*},u^{*},{\bf
p})\theta=-\int_{0}^{Z}\left[\theta(z,T)p(z,T)-\theta(z,0)p(z,0)\right]\,\mathrm{d}z\\\
\quad+\int_{Q}\theta\left[\frac{\partial p}{\partial
t}+(1-f(\theta^{*}))\frac{\mathrm{d}f}{\mathrm{d}\theta}(\theta^{*})+\beta\frac{\partial^{2}p}{\partial
z^{2}}+\frac{\mathrm{d}K}{\mathrm{d}\theta}(\theta^{*})\frac{\partial
p}{\partial z}\right]\,\mathrm{d}z\mathrm{d}t\\\
\quad+\int_{0}^{T}\left[\frac{\partial\theta}{\partial
z}(Z,t)\left(\beta(\theta^{*})p\right)_{\big{|}z=Z}-\frac{\partial\theta}{\partial
z}(0,t)\left(\beta(\theta^{*})p\right)_{\big{|}z=0}\right]\,\mathrm{d}t\\\
\quad+\int_{0}^{T}\theta(Z,t)\left[-\beta(\theta^{*})\frac{\partial
p}{\partial
z}-\frac{\mathrm{d}K}{\mathrm{d}\theta}(\theta^{*})p+\frac{\mathrm{d}\beta}{\mathrm{d}\theta}(\theta^{*})\frac{\partial\theta^{*}}{\partial
z}p-p_{2}\right]\,\mathrm{d}t\\\
\quad-\int_{0}^{T}\theta(0,t)\left[-\beta(\theta^{*})\frac{\partial
p}{\partial
z}-\frac{\mathrm{d}K}{\mathrm{d}\theta}(\theta^{*})p+\frac{\mathrm{d}\beta}{\mathrm{d}\theta}(\theta^{*})\frac{\partial\theta^{*}}{\partial
z}p+p_{1}\right]\,\mathrm{d}t\;.$

Thus, we deduce that the adjoint variable ${\bf p}$ satisfies

$\begin{cases}\frac{\partial p}{\partial
t}+\beta(\theta^{*})\frac{\partial^{2}p}{\partial
z^{2}}+\frac{\mathrm{d}K}{\mathrm{d}\theta}(\theta^{*})\frac{\partial
p}{\partial z}=F(\theta^{*})&\text{ in }Q,\\\ p(z,T)=0&\text{ in }\
(0,Z)\,,\\\ p(0,t)=p(Z,t)=0&\text{ for }t\in(0,T)\,,\\\
p_{1}=\left(\beta(\theta^{*})\frac{\partial p}{\partial
z}\right)_{\big{|}z=0}&\\\ p_{2}=-\left(\beta(\theta^{*})\frac{\partial
p}{\partial z}\right)_{\big{|}z=Z}&\end{cases}$

where

$F(\theta)\mathrel{\mathop{:}}=[f(\theta)-1]\frac{\mathrm{d}f}{\mathrm{d}\theta}(\theta).$

On the other hand, since

$\displaystyle D_{u}\mathcal{L}(\theta^{*},u^{*},{\bf p})u$
$\displaystyle=\int_{0}^{T}(\lambda u^{*}+p_{1})u\,\mathrm{d}t\;,$

the condition $D_{u}\mathcal{L}(\theta^{*},u^{*},p^{*})\cdot(u-u^{*})\geq 0$
for all $u\in U_{ad}$ implies the optimality condition

$\left\langle\lambda u^{*}(t)+\left(\beta(\theta^{*})\frac{\partial
p}{\partial z}\right)_{\big{|}z=0},{u}-u^{*}\right\rangle_{L^{2}(0,T)}\geq 0$

for all ${u}\in U_{ad}$. We thus obtain that any optimal solution
$(\theta^{*},u^{*},p^{*})$ of problem (2.1)-(4.1)-(4.2) must satisfy the
optimality system

(4.3a) $\displaystyle\begin{cases}\frac{\partial\theta}{\partial
t}-\frac{\partial^{2}\beta^{*}(\theta)}{\partial z^{2}}+\frac{\partial
K(\theta)}{\partial z}=f&\text{ in }Q,\\\ \theta(z,0)=\theta_{0}(z)&\text{ in
}(0,Z),\\\ \theta(0,t)=v(t)&\text{ for }t\in(0,T),\\\ \theta(Z,t)=g(t)&\text{
for }t\in(0,T),\\\ \end{cases}$ (4.3b)
$\displaystyle\begin{cases}\frac{\partial p}{\partial
t}+\beta(\theta^{*})\frac{\partial^{2}p}{\partial
z^{2}}+\frac{\mathrm{d}K}{\mathrm{d}\theta}(\theta^{*})\frac{\partial
p}{\partial z}=F(\theta^{*})&\text{ in }Q,\\\ p(z,T)=0&\text{ in }\
(0,Z)\,,\\\ p(0,t)=p(Z,t)=0&\text{ for }t\in(0,T)\,,\end{cases}$ (4.3c)
$\displaystyle\left\langle\lambda
u^{*}(t)+\left(\beta(\theta^{*})\frac{\partial p}{\partial
z}\right)_{\big{|}z=0},{u}-u^{*}\right\rangle_{L^{2}(0,T)}\geq 0$

for all ${u}\in U_{ad}$, where we recall that $v=\theta_{r}+u$. In the next
section, we exploit this optimality system to build suitable algorithms to
numerically solve the optimal control problem (4.1)-(2.1).

##  5\. Algorithm and numerical simulations

Our optimization procedure will follow the Projected Gradient Descent (PGD)
described in Algorithm 1 (see [35] for a thorough introduction to such
optimization algorithms). However, when solving (4.3) with PGD, it could
happen that Theorem 3.3 is not satisfied at each iteration, thus incurring
numerical difficulties due to the singularity of water diffusivity $\beta$ at
$\theta=\theta_{S}$. Therefore, we shall approximate it by truncation: given a
small $\varepsilon>0$, we define

(5.1)
$\beta_{\varepsilon}(r)\mathrel{\mathop{:}}=\begin{cases}\beta(r),&r\leq\theta_{S}-\varepsilon,\\\
\beta(\theta_{S}-\varepsilon),&r>\theta_{S}-\varepsilon,\end{cases}$

as shown in the Figure 1. Regularization (5.1) is a standard technique when
dealing with Richards’ equation to handle singularities in the diffusion term,
and it is used in finite difference schemes [30, 36] or FEM [37] for both the
mathematical and numerical analysis of degenerate, and possibly doubly-
degenerate, parabolic equations.
In the following simulations, we are then actually computing the numerical
solutions to (4.3) after replacing $\beta$ with $\beta_{\varepsilon}$ as
defined in (5.1).

Figure 1. Graph of the truncated water diffusivity $\beta_{\varepsilon}$ from
(5.1).

Moreover, we have set a maximum number of iterations equal to 100 before
exiting PGD iterations, a tolerance of $10^{-5}$, and a regularization
parameter $\varepsilon=10^{-3}$ for water diffusivity in (5.1). We stress that
the order of magnitude of $\varepsilon$ has been chosen so to be consistent
with that of the different $\theta_{S}$ values selected in all the simulations
that follow.

Algorithm 1 Projected Gradient Descent Algorithm applied to Richards’ equation
for determining optimal irrigation. Intermediate steps for solving direct
(line 3) and adjoint (line 4) problems and for computing the optimal descent
direction step-size (line 6) are performed using MATLAB pdepe and fmincon
functions, respectively.

1:$n\leftarrow 1,u_{n}\leftarrow u_{0}$$\triangleright$ Initialization

2:while $n\leq\textrm{maxit}$ do

3: $\theta_{n}\leftarrow$ directProblem($u_{n}$)$\triangleright$ Solving the
I-BV problem

4: $p_{n}\leftarrow$ adjointProblem($\theta_{n}$)$\triangleright$ Solving the
adjoint problem

5: $r_{n}\leftarrow-\left(\lambda
u_{n}+\beta(p_{n}(z=0))\partial_{z}p_{n}(z=0)\right)$$\triangleright$ Descent
direction

6:
$s_{n}\leftarrow\arg\min_{s}J(\mathrm{pr}_{\mathbf{U}}(u_{n}+sr_{n}))$$\triangleright$
Optimal direction Step-size

7:

8: if $\left\lvert
J(\mathrm{pr}_{\mathbf{U}}(u_{n}+s_{n}r_{n}))-J(u)\right\rvert<\mathrm{tol}$
then$\triangleright$ Exit condition

9: break

10: else

11: $u\leftarrow\mathrm{pr}_{\mathbf{U}}(u_{n}+s_{n}r_{n})$$\triangleright$
Updating control

12: $n\leftarrow n+1$

13: end if

14:end while

Moreover, we select a root water uptake model of Feddes type (used, for
instance, in [4, 6]) as source term in (2.1). Its expression is given by

(5.2)
$f(h)=\varphi\widehat{f}(h),\quad\widehat{f}(h)\mathrel{\mathop{:}}=\begin{cases}0,&\textrm{
if }h_{1}\leq h\leq 0\textrm{ or }h\leq h_{4},\\\
\frac{h-h_{1}}{h_{2}-h_{1}},&\textrm{if }h_{2}<h<h_{1},\\\
1,&\textrm{if\quad}h_{3}\leq h\leq h_{2},\\\
\frac{h-h_{4}}{h_{3}-h_{4}},&\textrm{if }h_{4}<h<h_{3},\end{cases}$

with the following values, in cm: $h_{4}\approx-820,\ h_{3}\approx-400,\
h_{2}\approx-350,\ h_{1}=0$. Also, we set $\varphi=0.1/Z$, where $Z$ is the
soil depth, and $\lambda=0.1$ in (4.1).

######  Remark 5.1.

Let us notice that the maximum value for $\widehat{f}(h)$ in (5.2) is set to
$1$ for normalization purposes. In fact, when it comes to practical problems,
one uses $f(h)$ properly rescaled according to experimental evidences through
the factor $\varphi$, which is the ratio of the potential transpiration rate
and the rooting depth, as explained in [38].
Moreover, we stress that in general one does not necessarily require the
source term $f(h)$ to be zero for the values of $h$ corresponding to the
boundary of $[\theta_{r},\theta_{S}]$. However, from a physical point of view,
it makes sense for a source term to vanish when the soil is either dry or
saturated. This is exactly the case of Feddes-type source terms as the one we
consider in (5.2) for our numerical simulations.

In Example 5.2 and Example 5.3 below, we simulate a soil described by
Haverkamp model [16, 39], whose constitutive relations are given by

(5.3)
$\theta(h)=\frac{\alpha\left(\theta_{S}-\theta_{r}\right)}{\alpha+\left\lvert
h\right\rvert^{\beta_{2}}}+\theta_{r}\;,\qquad
K(h)=K_{S}\frac{A}{A+\left\lvert h\right\rvert^{\beta_{1}}}\;,$

representing water retention curve and hydraulic conductivity, respectively.
We first verify that Haverkamp model falls within the quasi-unsaturated model
for a suitable choice of the parameters involved in the model setting. In
fact, in Haverkamp model we have

(5.4)
$\beta(\theta(h))=\frac{K_{S}A(\alpha+(-h)^{\beta_{2}})^{2}}{(A+(-h)^{\beta_{1}})\alpha(\theta_{S}-\theta_{r})\beta_{2}(-h)^{\beta_{2}-1}}$

with $\beta_{1},\beta_{2}>0$, that is Lipschitz, monotonically increasing and
bounded from below. In order to satisfy assumption ($\mathbf{H_{\beta}}$), we
shall ensure that

$\lim_{h\nearrow 0^{-}}\beta(\theta(h))=+\infty.$

From (5.4), a straightforward computation provides that this condition is
satisfied if and only if $\beta_{2}>1$.

######  Example 5.2.

This is the case of the sandy soil considered in [16], with parameters

(5.5) $\begin{split}K_{S}=34\,{\rm cm/h},\,\,A=1.175\times
10^{6},\,\,\beta_{1}=4.74,\\\
\theta_{S}=0.287,\,\,\theta_{r}=0.075,\alpha=1.611\times
10^{6},\,\,\beta_{2}=3.96,\end{split}$

where the root water uptake model is as in (5.2).

We have performed a simulation of $T=3$ hours with a maximum depth of $Z=70$
cm.

As in (2.1), boundary condition at the top varies in time according to the
irrigation strategy:

(5.6) $\theta(0,t)=\theta_{\rm top}(t)=u(t)+\theta_{r},$

while bottom condition has been chosen so to be constant over time:

(5.7) $\theta(Z,t)=\theta_{\rm
bottom}(t)=0.9\theta_{r}+0.1\theta_{S},\,\,t\in[0,T].$

Finally, initial condition is linearly varying over time as

(5.8) $\theta(z,0)=\theta_{\rm top}(0)+z\frac{\theta_{\rm
bottom}(0)-\theta_{\rm top}(0)}{Z},\quad z\in[0,Z].$

(a)

(b)

(c)

Figure 2. Numerical simulations relative to Example 5.2, where initial
condition is given by (5.8) and boundary conditions are as in (5.6) and (5.7),
respectively.

Our simulations have been produced using MATLAB active-set algorithm; we
report that same results are obtained using MATLAB sqp. We have observed that
convergence is reached after 3 iterates within the given tolerance, and the
numerical solution is locally optimal. Results are in Figure 2. As can be
seen, the optimal control framework succeeds in determining an optimal control
that optimizes the performance index (4.1), with a reduced water consumption
and average water content over time.

######  Example 5.3.

In this second simulation, using the same soil as in Example 5.2, we consider
a time-varying bottom condition

(5.9)
$\theta(Z,t)=\theta_{\textrm{bottom}}(t)=\left(1-\frac{t}{T}\right)\theta_{b1}+\frac{t}{T}\theta_{b2},\,\,t\in[0,T],$

where

$\theta_{b1}\mathrel{\mathop{:}}=0.9\theta_{r}+0.1\theta_{S},\,\theta_{b2}\mathrel{\mathop{:}}=0.7\theta_{r}+0.3\theta_{S},$

while top condition is given by

(5.10) $\theta(0,t)=\theta_{\textrm{top}}(t)=u(t)+\theta_{r}.$

Moreover, initial condition is

(5.11) $\theta(z,0)=\theta_{\rm top}(0)+z\frac{\theta_{\rm top}(0)-\theta_{\rm
bottom}(0)}{Z},\quad z\in[0,Z].$

It turns out that MATLAB sqp converges, in $3$ iterates, to a local optimal
solution, further providing the best results if compared to active-set.
Results are displayed in Figure 3.

(a)

(b)

(c)

Figure 3. Numerical simulations relative to Example 5.3. For this simulation,
initial condition is given by (5.11) and boundary conditions are given by
(5.9) and (5.10), respectively.

In Example 5.4 and Example 5.5 that follow, we consider the classical Van
Genuchten-Mualem constitutive relations in the unsaturated zone, given by

(5.12a) $\displaystyle\theta\left(\psi\right)$
$\displaystyle=\theta_{r}+\frac{\theta_{S}-\theta_{r}}{\left(1+\left\lvert\alpha\psi\right\rvert^{n}\right)^{m}},\quad
m\mathrel{\mathop{:}}=1-\frac{1}{n},$ (5.12b) $\displaystyle k(\psi)$
$\displaystyle=K_{S}\left[\frac{1}{1+\left\lvert\alpha\psi\right\rvert^{n}}\right]^{\frac{m}{2}}\left[1-\left(1-\frac{1}{1+\left\lvert\alpha\psi\right\rvert^{n}}\right)^{m}\right]^{2}.$

In order to verify under which conditions Van Genuchten-Mualem model satisfies
the quasi-unsaturated model, we need to analyze its corresponding function
$\beta(\theta(h))$. Letting

$\varphi(h)\mathrel{\mathop{:}}=\frac{1}{1+\left\lvert\alpha
h\right\rvert^{n}},$

from (5.12) and exploiting the fact that $h<0$, it follows that

(5.13)
$\beta(\theta(h))=\frac{K_{S}\left[1-(1-\varphi(h))^{m}\right]^{2}}{mn\alpha^{n}(\theta_{S}-\theta_{r})\left\lvert
h\right\rvert^{n-1}\varphi(h)^{\frac{m}{2}+1}}.$

It is an easy computation that $\beta$ is Lipschitz, monotonically increasing
and bounded from below. In order to satisfy assumption ($\mathbf{H_{\beta}}$),
there needs

$\lim_{h\nearrow 0^{-}}\beta(\theta(h))=+\infty.$

From (5.13), this condition is satisfied if and only if $n>1$.

This is the case for the simulations reported below. More specifically, we are
going to consider a Berino loamy fine sand and a Glendale clay loam, with
parameters drawn from [40, 26].

######  Example 5.4.

The Berino loamy fine sand is defined by the following hydraulic parameters:

$\theta_{r}=0.0286,\theta_{S}=0.3658,\alpha=0.0280,n=2.2390,K_{S}=22.5416\,\
\mathrm{cm/h}.$

Here, as in Example 5.3, we consider a time-varying bottom condition

(5.14)
$\theta(Z,t)=\theta_{\textrm{bottom}}(t)=\left(1-\frac{t}{T}\right)\theta_{b1}+\frac{t}{T}\theta_{b2},\,\,t\in[0,T],$

where

$\theta_{b1}\mathrel{\mathop{:}}=0.3\theta_{r}+0.7\theta_{S},\,\theta_{b2}\mathrel{\mathop{:}}=0.1\theta_{r}+0.9\theta_{S};$

boundary condition at the top of the domain is, again as in Example 5.3,

(5.15) $\theta(0,t)=\theta_{\textrm{top}}(t)=u(t)+\theta_{r}.$

However, now initial condition is a quadratic polynomial function of depth and
is set as

(5.16) $\theta(z,0)\mathrel{\mathop{:}}=(\theta_{\rm bottom}(0)-\theta_{\rm
top}(0))\left(\frac{z}{Z}\right)^{2}+\theta_{\rm top}(0),\quad z\in[0,Z].$

Results relative to this soil are depicted in Figure 4 and are obtained using
MATLAB active-set, where $Z=50$ cm and $T=12$ hours.

######  Example 5.5.

Simulations on the Glendale clay loam are obtained using the following
parameters:

$\theta_{r}=0.1060,\theta_{S}=0.4686,\alpha=0.0104,n=1.3954,K_{S}=0.5458\,\,\mathrm{cm/h}.$

For this experiment, we fix $Z=30$ cm and $T=36$ hours; bottom boundary
condition is given by

(5.17)
$\theta(Z,t)=\theta_{\textrm{bottom}}(t)=\left(1-\frac{t}{T}\right)\theta_{b1}+\frac{t}{T}\theta_{b2},\,\,t\in[0,T],$

where

$\theta_{b1}\mathrel{\mathop{:}}=0.5\theta_{r}+0.5\theta_{S},\,\theta_{b2}\mathrel{\mathop{:}}=0.7\theta_{r}+0.3\theta_{S},$

and top boundary condition is, as in previous examples,

(5.18) $\theta(0,t)=\theta_{\textrm{top}}(t)=u(t)+\theta_{r}.$

Initial condition is again set as

(5.19) $\theta(z,0)\mathrel{\mathop{:}}=(\theta_{\rm bottom}(0)-\theta_{\rm
top}(0))\left(\frac{z}{Z}\right)^{2}+\theta_{\rm top}(0),\quad z\in[0,Z].$

Results are depicted in Figure 5. Here, we employed to MATLAB sqp for solving
the optimization problem by PGD.

(a)

(b)

(c)

Figure 4. Numerical simulations relative to Berino loamy fine sand in Example
5.4, where we set bottom and top boundary conditions as in (5.14) and (5.15),
respectively, whilst cubic initial condition is as in (5.16).

(a)

(b)

(c)

Figure 5. Numerical simulations relative to Glendale clay loam in Example 5.5.
Bottom condition is given by (5.17), top condition by (5.18); initial
condition is a quadratic polynomial as defined in (5.19).

##  6\. Conclusions

In this paper we introduce an optimal control approach aimed at optimizing the
water content provided by irrigation, applying Richards’ equation for
unsaturated flow. We make use of quasi-unsaturated model introduced in [33],
extending the well-posedness results for nonlinear sink terms and deriving
suitable optimality conditions for an irrigation performance index of tracking
type.
We set the model within a MATLAB solver by implementing a properly adapted
Projected Gradient Descent method, and provide significant numerical results
over a meaningful variety of soils; a deeper analytical treatise of the
control system is beyond the scopes of this paper, and it is currently under
investigations by the authors.
This paper could pave the way to an extensive use of control techniques for
optimizing irrigation in real life applications, and this framework could
easily be incorporated in existing irrigation software based on Richards’
equation solvers. Moreover, it is worth investigating qualitative features of
the more general _saturated-unsaturated_ model, for which there is an
increasing need of both numerical and analytical results and approaches. In
this context, tools from set-valued analysis and discrete control techniques
could carry improvements in understanding such problems.

## Acknowledgments

MB acknowledges the partial support of RIUBSAL project funded by Regione
Puglia under the call “P.S.R. Puglia 2014/2020 - Misura 16 – Cooperazione -
Sottomisura 16.2 “Sostegno a progetti pilota e allo sviluppo di nuovi
prodotti, pratiche, processi e tecnologie”: in particular he thanks Mr.
Giuseppe Leone and Mrs. Gina Dell’Olio for supporting the project activities;
FVD has been supported by REFIN Project, grant number 812E4967, funded by
Regione Puglia: both authors acknowledge the partial support of GNCS-INdAM. RG
acknowledges the support of the Natural Sciences and Engineering Research
Council of Canada (NSERC), funding reference number RGPIN-2021-02632.

## References

  * [1] Y. Mao, S. Liu, J. Nahar, J. Liu, F. Ding, Soil moisture regulation of agro-hydrological systems using zone model predictive control, Computers and Electronics in Agriculture 154 (2018) 239–247.  doi:https://doi.org/10.1016/j.compag.2018.09.011.
  * [2] A. Coppola, G. Dragonetti, A. Sengouga, N. Lamaddalena, A. Comegna, A. Basile, N. Noviello, L. Nardella, Identifying optimal irrigation water needs at district scale by using a physically based agro-hydrological model, Water 11 (4) (2019).  doi:10.3390/w11040841.
  * [3] W. R. Gardner, Modeling water uptake by roots, Irrigation Science 12 (3) (1991) 109–114.  doi:10.1007/BF00192281.
  * [4] F. V. Difonzo, C. Masciopinto, M. Vurro, M. Berardi, Shooting the numerical solution of moisture flow equation with root uptake: a Python tool, Water Resources Management 35 (2021) 2553–2567.  doi:10.1007/s11269-021-02850-2.
  * [5] E. Nazari, S. Besharat, K. Zeinalzadeh, A. Mohammadi, Measurement and simulation of the water flow and root uptake in soil under subsurface drip irrigation of apple tree, Agricultural Water Management 255 (2021) 106972.  doi:https://doi.org/10.1016/j.agwat.2021.106972.
  * [6] S. K. Deb, M. K. Shukla, J. Šimůnek, J. G. Mexal, Evaluation of Spatial and Temporal Root Water Uptake Patterns of a Flood-Irrigated Pecan Tree Using the HYDRUS (2D/3D) Model, Journal of Irrigation and Drainage Engineering 139 (8) (2013) 599–611.  doi:10.1061/(ASCE)IR.1943-4774.0000611.
  * [7] S. O. Lopes, F. A. C. C. Fontes, R. M. S. Pereira, M. de Pinho, A. M. Gonçalves, Optimal control applied to an irrigation planning problem, Mathematical Problems in Engineering 2016 (2016) 5076879.  doi:https://doi.org/10.1155/2016/5076879.
  * [8] M. Berardi, M. D’Abbicco, G. Girardi, M. Vurro, Optimizing water consumption in Richards’ equation framework with step-wise root water uptake: a simplified model, Transport in Porous Media 142 (2022) 469–498.  doi:https://doi.org/10.1007/s11242-021-01730-y.
  * [9] N. I. Challapa Molina, J. P. V.S. Cunha, Non-collocated sliding mode control of partial differential equations for soil irrigation, Journal of Process Control 73 (2019) 1–8.  doi:https://doi.org/10.1016/j.jprocont.2018.11.002.
  * [10] F. Wein, N. Chen, N. Iqbal, M. Stingl, M. Avila, Topology optimization of unsaturated flows in multi-material porous media: Application to a simple diaper model, Communications in Nonlinear Science and Numerical Simulation 78 (2019) 104871.  doi:https://doi.org/10.1016/j.cnsns.2019.104871.
  * [11] L. Bergamaschi, M. Putti, Mixed finite elements and Newton-type linearizations for the solution of Richards’equation, International Journal for Numerical methods in Engineering 45 (1999) 1025–1046.  doi:10.1002/(SICI)1097-0207(19990720)45:8<1025::AID-NME615>3.0.CO;2-G.
  * [12] V. Casulli, P. Zanolli, A Nested Newton-Type Algorithm for Finite Volume Methods Solving Richards’ Equation in Mixed Form, SIAM Journal on Scientific Computing 32 (4) (2010) 2255–2273.  doi:10.1137/100786320.
  * [13] I. Pop, F. Radu, P. Knabner, Mixed finite elements for the Richards’ equation: linearization procedure, Journal of Computational and Applied Mathematics 168 (1–2) (2004) 365 – 373.  doi:https://doi.org/10.1016/j.cam.2003.04.008.
  * [14] F. List, F. A. Radu, A study on iterative methods for solving Richards’ equation, Computational Geosciences 20 (2) (2016) 341–353.  doi:10.1007/s10596-016-9566-3.
  * [15] K. Mitra, I. Pop, A modified L-scheme to solve nonlinear diffusion problems, Computers & Mathematics with Applications (2018).  doi:https://doi.org/10.1016/j.camwa.2018.09.042.
  * [16] M. A. Celia, E. T. Bouloutas, R. L. Zarba, A general mass-conservative numerical solution for the unsaturated flow equation, Water Resources Research 26 (7) (1990) 1483–1496.  doi:10.1029/WR026i007p01483.
  * [17] T. Arbogast, M. Wheeler, N. Zhang, A Nonlinear Mixed Finite Element Method for a Degenerate Parabolic Equation Arising in Flow in Porous Media, SIAM Journal on Numerical Analysis 33 (4) (1996) 1669––1687.  doi:10.1137/S0036142994266728.
  * [18] E. Schneid, P. Knabner, F. Radu , A priori error estimates for a mixed finite element discretization of the Richards’ equation, Numerische Mathematik 98 (2004) 353–370.  doi:https://doi.org/10.1007/s00211-003-0509-2.
  * [19] C. Kees, M. Farthing, C. Dawson, Locally conservative, stabilized finite element methods for variably saturated flow, Computer Methods in Applied Mechanics and Engineering 197 (51) (2008) 4610–4625.  doi:https://doi.org/10.1016/j.cma.2008.06.005.
  * [20] H. Li, M. Farthing, C. Miller, Adaptive local discontinuous Galerkin approximation to Richards’ equation, Advances in Water Resources 30 (9) (2007) 1883–1901.  doi:https://doi.org/10.1016/j.advwatres.2007.02.007.
  * [21] J.-B. Clément, F. Golay, M. Ersoy, D. Sous, An adaptive strategy for discontinuous Galerkin simulations of Richards’ equation: Application to multi-materials dam wetting, Advances in Water Resources 151 (2021) 103897\.  doi:https://doi.org/10.1016/j.advwatres.2021.103897.
  * [22] R. Eymard, M. Gutnic, D. Hilhorst, The finite volume method for Richards equation, Computational Geosciences 3 (3-4) (1999) 259–294.  doi:10.1023/A:1011547513583.
  * [23] G. Manzini, S. Ferraris, Mass-conservative finite volume methods on 2-d unstructured grids for the Richards’ equation, Advances in Water Resources 27 (12) (2004) 1199 – 1215.  doi:https://doi.org/10.1016/j.advwatres.2004.08.008.
  * [24] D. Seus, K. Mitra, I. S. Pop, F. A. Radu, C. Rohde, A linear domain decomposition method for partially saturated flow in porous media, Computer Methods in Applied Mechanics and Engineering 333 (2018) 331 – 355.  doi:https://doi.org/10.1016/j.cma.2018.01.029.
  * [25] T. Hoang, I. Pop, Iterative methods with nonconforming time grids for nonlinear flow problems in porous media, Acta Mathematica Vietnamica (2022).  doi:https://doi.org/10.1007/s40306-022-00486-x.
  * [26] M. Berardi, F. Difonzo, M. Vurro, L. Lopez, The 1D Richards’ equation in two layered soils: a Filippov approach to treat discontinuities, Advances in Water Resources 115 (2018) 264–272.  doi:10.1016/j.advwatres.2017.09.027.
  * [27] M. Berardi, F. V. Difonzo, L. Lopez, A mixed MoL-TMoL for the numerical solution of the 2D Richards’ equation in layered soils, Computers & Mathematics with Applications 79 (2020) 1990–2001.  doi:https://doi.org/10.1016/j.camwa.2019.07.026.
  * [28] M. W. Farthing, F. L. Ogden, Numerical solution of Richards’ equation: A review of advances and challenges, Soil Science Society of America Journal 81 (8) (2017) 04017025.  doi:doi:10.2136/sssaj2017.02.0058.
  * [29] Y. Zha, J. Yang, J. Zeng, C.-H. M. Tso, W. Zeng, L. Shi, Review of numerical solution of Richardson–Richards equation for variably saturated flow in soils, WIREs Water 6 (5) e1364.  doi:https://doi.org/10.1002/wat2.1364.
  * [30] H. W. Alt, S. Luckhaus, Quasilinear elliptic-parabolic differential equations, Mathematische Zeitschrift 183 (3) (1983) 311–341.  doi:10.1007/BF01176474.
  * [31] F. Otto, L1-Contraction and Uniqueness for Quasilinear Elliptic–Parabolic Equations, Journal of Differential Equations 131 (1) (1996) 20–38.  doi:https://doi.org/10.1006/jdeq.1996.0155.
  * [32] B. Schweizer, Regularization of outflow problems in unsaturated porous media with dry regions, Journal of Differential Equations 237 (2) (2007) 278–306.  doi:https://doi.org/10.1016/j.jde.2007.03.011.
  * [33] G. Marinoschi, Functional Approach to Nonlinear Models of Water Flow in Soils, Springer, Dordrecht, The Netherlands, 2006.
  * [34] W. Merz, P. Rybka, Strong solutions to the Richards equation in the unsaturated zone, Journal of Mathematical Analysis and Applications 371 (2) (2010) 741–749.  doi:https://doi.org/10.1016/j.jmaa.2010.05.066.
  * [35] F. Tröltzsch, J. Sprekels, Optimal Control of Partial Differential Equations: Theory, Methods, and Applications, Graduate studies in mathematics, American Mathematical Society, Providence R.I, 2010.
  * [36] I. S. Pop, B. Schweizer, Regularization schemes for degenerate richards equations and outflow conditions, Mathematical Models and Methods in Applied Sciences 21 (2011) 1685–1712.
  * [37] R. H. Nochetto, A. Schmidt, C. Verdi, Adapting meshes and time-steps for phase change problems, Atti della Accademia Nazionale dei Lincei. Classe di Scienze Fisiche, Matematiche e Naturali. Rendiconti Lincei. Matematica e Applicazioni 8 (4) (1997) 273–292.
  * [38] A. Utset, M. E. Ruiz, J. Garcia, R. A. Feddes, A swacrop-based potato root water-uptake function as determined under tropical conditions, Potato Research 43 (1) (2000) 19–29.  doi:10.1007/BF02358510.
  * [39] M. Berardi, F. Difonzo, F. Notarnicola, M. Vurro, A transversal method of lines for the numerical modeling of vertical infiltration into the vadose zone, Applied Numerical Mathematics 135 (2019) 264–275.  doi:https://doi.org/10.1016/j.apnum.2018.08.013.
  * [40] R. G. Hills, I. Porro, D. B. Hudson, P. J. Wierenga, Modeling one-dimensional infiltration into very dry soils: 1. model development and evaluation, Water Resources Research 25 (6) (1989) 1259–1269.  doi:10.1029/WR025i006p01259.


Heliosphysics 2024 Science White Paper Heliosphere Meets Interstellar Medium,
in a Galactic Context

Category: Basic Research
Primary Topic: Emerging Opportunities (Interstellar Medium and Astrospheres)
Secondary Topic: Outer Heliosphere

Principal Author:

Name: Stella Koch Ocker Institution: Cornell University Email:
<EMAIL_ADDRESS>
Co-authors: James Cordes(1), Shami Chatterjee(1), Jeffrey Hazboun(2,3),
Timothy Dolch(4,5), Daniel Stinebring(6), Dustin Madison(7), Stephen White(8),
Gregory Taylor(9), Natalia Lewandowska(10), Michael Lam(11) (1)Cornell
University, (2)Oregon State University, (3)University of Washington Bothell,
(4)Hillsdale College, (5)Eureka Scientific Inc., (6)Oberlin College,
(7)University of the Pacific, (8)Air Force Research Lab, (9)University of New
Mexico, (10)State University of New York Oswego, (11)Rochester Institute of
Technology

Synopsis: The physical conditions within our heliosphere are driven by the
Sun’s motion through an evolving interstellar environment that remains largely
unexplored. The next generation of outer heliosphere and interstellar
explorers will answer fundamental questions about the heliosphere’s
relationship with the very local interstellar medium (VLISM) by diving deeper
into the Sun’s interstellar surroundings. The impact of these future missions
will be vastly enhanced by concurrent, interdisciplinary studies that examine
the direct connections between conditions within the heliosphere, the
heliosphere’s immediate interstellar environment, and the larger-scale
Galactic ISM. Comparisons of the heliosphere and VLISM to their analogs across
the Galaxy will constrain the global processes shaping both stellar
astrospheres and their sustained impact on the ISM.

Figure 1: a) Galactic dust polarization map [1]. b) The Galaxy in H$\alpha$
[2]. c) Bow shock of the star LL Ori (credit: NASA/ESA). d) Last image of
Earth from Voyager 1 (credit: NASA/JPL).

##  1 Introduction

The mid-2012 arrival of Voyager 1 in the VLISM [3], the region of interstellar
space just outside of the heliosphere, constituted the beginning of a new era
of heliophysics research. In just 10 years, Voyager 1 (V1), and later Voyager
2 (V2), have demonstrated that the heliopause is a highly non-uniform, porous
boundary between the heliosphere and VLISM [4, 5, 6, 7]. Independent
astronomical observations of the local interstellar cloud (LIC) and the Local
Bubble (LB) [8, 9], combined with magnetohydrodynamic (MHD) modeling of the
heliosphere [10], have begun to contextualize Voyager’s findings in terms of
the long-term journey and evolution of the heliosphere’s global shape, size,
and energetics. Despite this progress, critical questions about the
heliosphere’s relationship with the ISM remain unanswered, and even
exacerbated. In the next decade, a combination of outer heliosphere and VLISM
missions with astronomical observations of the local ISM and other stars will
answer fundamental, interconnected questions:

  * •

What is the nature and extent of the Sun’s interaction with the ISM, and how
does this relationship compare to other stars?

  * •

Can a spacecraft explorer travel beyond the VLISM, and if so, what will it
reveal about the content and energetics of the Galactic ISM?

  * •

How unusual or typical is the Sun’s interstellar environment, compared to that
of other astrospheres throughout the Galaxy?

Addressing these questions requires the assessment of both stochastic
(turbulent) and deterministic processes near stellar-ISM interfaces, in
comparison to analogous processes in the diffuse, larger-scale ISM. Novel
interdisciplinary approaches that connect spacecraft data to remote
observations of other stars and the ISM will illuminate not only the processes
driving the global structure of the heliosphere, but also the extent to which
these processes are common features of stellar systems across the Galaxy.

We focus on two key themes: spacecraft exploration of the VLISM, and
complementary studies of the heliosphere, local ISM, and other star-ISM
interactions enabled by observations of high-velocity stars and compact radio
sources, including pulsars and extragalactic quasars that act as beacons
through all ionized components of the ISM.

##  2 The Dynamic Heliosphere-ISM Interaction

Necessity of Interstellar Spacecraft. The Voyager Interstellar Mission (VIM)
is currently the only source of direct measurements in the VLISM. These in
situ data are incomparable to any other indirect probe of the outer
heliosphere and VLISM, and have led to key discoveries about their
relationship, including:

  * •

The magnetic boundary of the heliopause is not uniform. While both V1 and V2
traversed comparable plasma density discontinuities at the heliopause [11], V1
saw dramatic magnetic field fluctuations [4] that were not later observed by
V2 [6]. V2 also observed a cosmic ray boundary layer outside the heliopause,
which was not seen by V1 [7]. Immediately outside the heliopause, V2 found
VLISM plasma about $2\times$ hotter than expected ($\sim 40,000$ K) [12].

  * •

Shocks of solar origin propagate beyond the heliopause, generating plasma
oscillation events, magnetic field jumps, and cosmic ray bursts [13]. While
these shocks were observed regularly (about once a year) by V1 from the
heliopause at 121 au out to 145 au [13], no such events have been observed
since. It is unclear whether the recent absence of shocks is related to the
solar cycle or whether it signifies the extent of solar influence on the
VLISM.

  * •

The unexpected detection of faint, persistent plasma wave emission by V1 has
enabled a highly sampled density map of the VLISM [14], and demonstrates that
the VLISM electron distribution has a supra-thermal tail [15, 16]; see Figure
2.

  * •

The V1 observation of a large magnetic pressure front and density jump at
about 150 au [17] without any accompanying plasma oscillations raises the
question of whether V1 has identified another boundary layer in the VLISM.

Figure 2: a) VLISM electric field intensity in frequency (kHz) vs. time
(years) measured by the Voyager 1 Plasma Wave Subsystem. Bright patches of
power are solar shock-triggered plasma oscillation events. b) Same as a) but
smoothed in time and shown with a restricted color stretch. The narrow line
between the solar-origin plasma oscillation events corresponds to the
persistent plasma wave emission that enables a near-continuous sampling of the
plasma density. This figure was adapted from Ref.[14]. Interstellar Probe will
be capable of resolving this narrow plasma line with much higher signal-to-
noise and over a greater extent of the VLISM.

These discoveries have raised questions with significant implications for
conditions within the heliosphere: What drives the non-uniformity of the outer
heliosphere’s structure? How far does the VLISM extend? Does the heliosphere
have a bow shock (and if so, what is its nature)? Do shock waves from other
stars perturb the VLISM?

Direct spacecraft measurements are critical to answering these questions.
While VIM is expected to continue operating for the next few years, the launch
of a new Interstellar Probe mission [18, 19] in the next decade would be
pivotal for heliophysics and VLISM research.

Interstellar Probe will provide the first high-resolution map of plasma
densities, magnetic fields, and energetic particles spanning both the inner
and outer heliosphere and the VLISM (out to at least 300 au) [18, 19].
Anticipated results include: The effect of solar cycle on the size and
structure of the heliosphere and its boundary with the VLISM; signatures of
interstellar-origin shocks and structures; the evolution of turbulence and its
degree of anisotropy from outer heliosphere to VLISM; and constraints on the
heliospheric bow shock’s existence. By exploring a novel trajectory through
the heliosphere, Interstellar Probe will empirically test MHD models for the
global structure of the heliosphere and its evolution as a function of
interstellar environment.111See the Interstellar Probe white paper (lead
author Pontus Brandt) for more detail about the mission concept.

The scientific impact of Interstellar Probe will extend beyond heliophysics to
interstellar astrophysics research. Spacecraft measurements of magnetic field
and electron density fluctuations are the only direct probes of interstellar
turbulence, which is characterized using magnetic field and plasma density
power spectra (see Figure 3). V1 has thus far provided the only constraints on
the interstellar turbulence power spectrum at kinetic and dissipative scales
($\sim$ kilometers), which cannot be accessed with other remote probes of the
ISM (Figure 3). However, there is a significant gap in the turbulence spectrum
probed by VIM between length scales of $\sim 10^{3}$ km and $\sim 0.01$ au.
Interstellar Probe will span this gap, thereby providing the first self-
consistently measured turbulence spectrum across 10 orders of magnitude in
length scale.

Figure 3: The lefthand plot shows the composite power spectrum of turbulent,
interstellar electron density fluctuations as a function of length scale (au),
inferred from three different in situ and path-integrated probes illustrated
on the right: A. In situ H$\alpha$ images of the pulsar bow shock known as the
Guitar Nebula; B. Plasma densities measured in situ using plasma oscillations
detected by the Voyager 1 Plasma Wave Subsystem (PWS); C. Path-integrated
pulsar dispersion measure (DM) variations measured by the NANOGrav
Collaboration over months to years. The proposed Interstellar Probe will
sample the as-yet unexplored region indicated by the red question mark.

The Unique Value of Astronomical Remote Sensing. Despite their extraordinary
resolving power, spacecraft measurements are limited by the singular
sightlines and distances they can probe. Remote observations along many lines-
of-sight (LOSs) across the sky are thus required to contextualize in situ
spacecraft measurements, both in and outside the heliosphere.

Remote sensing observations can be categorized as either path-integrated or in
situ. Here we focus mainly on radio pulsars222Optical emission and absorption
lines are also critical path-integrated probes of the VLISM and local clouds;
see the white paper led by Linsky and Redfield., which undergo dispersion in
cold plasma. The resulting frequency-dependent time delay in their pulse
times-of-arrival depends on the path-integrated electron column density, i.e.
the dispersion measure (DM). Another quantity of interest is the Faraday
rotation measure (RM) from polarization measurements, which gives information
about the path-integrated magnetic field. Together, DMs and RMs for many LOSs
are inputs to tomographic maps of the ionized ISM. Pulsars and other compact
radio sources are also scattered by intervening plasma, which is another path-
integrated effect discussed further below. The chief remote sensing method for
in situ observation is imaging. Imaging other stellar bow shocks and
astrospheres provides in situ constraints on interstellar gas and star-ISM
interactions that are complementary to spacecraft measurements of the
heliosphere and VLISM (see Section 3).

Inside the Heliosphere: Using Pulsar Dispersion through the Solar Wind.
Regular monitoring of a relatively dense pulsar population across the sky will
provide 3D constraints on the solar wind density throughout the heliosphere,
in a manner independent and complementary to the myriad space missions
studying solar weather. Pulsar DMs undergo annual variations as their LOSs
pass through the solar corona over the course of Earth’s orbit. These DM
variations are sensitive to solar wind features that are more complex than the
spherically symmetric $1/r^{2}$ component [20, 21]. Recent work [22, 23, 24,
25] has shown that pulsar timing arrays (PTAs) can be valuable probes of the
solar wind and its temporal behavior. PTAs survey a much wider swath of the
Sun’s environment than spacecraft by taking electron column density
information for 70+ LOSs every month, thereby probing even the outermost
reaches of the solar environment. Lower frequency, high-cadence observing
campaigns with, e.g., the Canadian Hydrogen Intensity Mapping Experiment
(CHIME) [26], the Low Frequency Array (LOFAR) [27, 28], and the Long
Wavelength Array (LWA) [[]]lwa_pulsar2015 are also continually adding data
that is extremely useful for DM variability studies, since these lower
frequencies allow for more accurate DM variation measurements [30]. Similar,
high-cadence monitoring of RMs will provide complementary magnetic field
information [31].

All of these data can add important tests to the well developed solar weather
modeling efforts undertaken by the heliophysics community [32]. These large-
scale measurements would allow for long-timescale monitoring of the fast and
slow solar wind, giving access to continuous measurements of the solar wind at
all solar latitudes. Solar wind models are also a critical element of high-
accuracy pulsar astronomy applications, including nanohertz frequency
gravitational wave detection (see Astro2020 decadal survey priority area “New
Messengers and New Physics”)333See also the Nanohertz Observatory for
Gravitational Waves (NANOGrav) Astro2020 White Paper:
https://baas.aas.org/pub/2020n7i195/release/1 [25].
Outside the Heliosphere: Interstellar Scintillation as Probes of Sub-au
Structure in the Local Interstellar Clouds (LICs). Radio scintillation
observations are uniquely sensitive to small-scale structure in the LICs that
shape the heliosphere. Radio emission from compact sources like pulsars and
extragalactic quasars is scattered by sub-au scale plasma inhomogeneities
along their LOSs, leading to constructive and destructive interference between
the observed scattered rays. The resulting diffraction pattern changes over
time due to the Earth’s motion relative to the LOS, a phenomenon known as
scintillation. Fourier analyses of pulsar scintillation patterns, a method
known as scintillometry, show that scattering appears to occur in confined,
discrete regions of plasma along the LOS (see Figure 4), and the precise
distances of these scattering regions are constrained by scintillation
observations taken months to years apart [33]. The scintillometry method
traces processes in ISM plasma from scales of $\sim 1000$ km to $\sim$au, thus
representing one of the few remote, path-integrated probes of sub-au scale
processes that are otherwise measured by spacecraft.

A number of pulsars exhibit scintillation from plasma structures within or
near the edge of the Local Bubble [34], and the presence of multiple
scattering regions along some pulsar LOSs suggests these structures have a
high number density in the ISM. In one case, the closest scattering region
lies just 5 pc away from Earth [33], within the LICs whose pressure,
temperature, density, and magnetic field determine the heliosphere’s current
shape and size and the flow of energetic particles past Earth [9] (see Figure
4).

Large-scale scintillometry of hundreds of pulsars with new, more capable
telescopes (e.g., DSA-2000444DSA-2000 Astro2020 White Paper:
https://baas.aas.org/pub/2020n7i255/release/1) will reveal more of these
structures, yielding unique information about the formation and sustenance of
micro-structure in the local interstellar clouds, including boundary
conditions on cloud interfaces, which are poorly understood. Dedicated pulsar
timing programs (with CHIME, MeerKAT, DSA-2000, and more) targeting large
numbers of nearby pulsars are required to achieve this goal, and are already
demonstrated to be highly successful at detecting scintillation from a large
fraction of such pulsars when observed [35, 34].

[capbesideposition=right,center,capbesidewidth=5.5cm]figure[]

Figure 4: a) Pulsar radio intensity in time vs. frequency, i.e. the dynamic
spectrum (from Ref.[36]). b) The scintillation arc spectrum formed via Fourier
transform of the dynamic spectrum. Each parabolic arc corresponds to a
distinct scattering screen. c) Locations of scattering screens along the LOS
to this pulsar, inferred with scintillometry (adapted from Ref.[33]). The
closest screen is 5 pc from Earth.

From the Heliosphere to the Local Bubble and ISM. The power of pulsar
measurements lies in precision DMs combined with precision parallaxes, which
yield an unbiased tracer of the ISM electron column density. Interpretation of
DMs does not require knowledge of the gas temperature or species abundances.
Telescopes at the lowest frequencies (e.g. the LWA) will produce the highest-
precision DMs possible for many pulsars. The current sample of pulsar
parallaxes within 400 pc of the heliosphere suggests the mean electron density
within the LB ranges from 0.009 to 0.018 cm-3 (Figure 5), calling into
question the pervasiveness of hot gas within the Bubble [37].

Hundreds of new pulsar parallaxes (from pulsar timing and very long baseline
interferometry) coupled with DMs will yield a high-resolution, tomographic
electron density map of the LB and immediate ISM, which is informative about
the ionizing radiation field that sustains the LICs. The LB represents just
one of many supernova-driven cavities that permeate the Galactic disk. Over
the next decade, up to 5,000 pulsars will be discovered, enabling
identification and characterization of the overall porosity in the ISM.

[capbesideposition=right,center,capbesidewidth=5cm]figure[]

Figure 5: Lines of sight to pulsars within 500 pc of the Sun with precisely
measured parallax distances, mapped on to the Galactic plane (the Galactic
Center is to the right). Boxes next to each point list the mean electron
density using the DM and parallax distance. Red lines indicate the three
objects closer than 200 parsecs.

##  3 The Galactic Context

Lessons from Stellar Bow Shocks. Observations of high-velocity stellar bow
shocks present a critical opportunity to constrain in real-time the ISM
conditions that also affect the astrospheres of lower velocity stars,
including the heliosphere. Bow shocks are observed around stars at a range of
life stages, including runaway OB stars and supergiants, and their signatures
span radio to X-ray wavelengths [38, 39, 40]. Neutron stars provide one test
bed for examining the effects of ISM conditions on bow shock evolution, as
they are born at speeds ranging from 100s to 1000s of kilometers per second
and generally produce bow shocks once they have moved from their supernova
remnants into the ISM. One of the few direct methods for detecting neutron
star bow shocks is by observing the H$\alpha$ emission produced by collisional
excitation of interstellar gas at the bow shock. Figure 3 shows Hubble Space
Telescope images (H$\alpha$) of the pulsar bow shock known as the Guitar
Nebula [41]. The bow shock morphology evolves over years-long timescales as a
direct result of ISM density variations sampled by the rapidly moving pulsar
[42, 43].

Similar ISM variations in plasma density, ionization, and magnetic fields
directly impact the flux of cosmic rays and pick-up ions on Earth and the
other planets [44]. As the heliosphere travels through a complex system of
interstellar clouds with varying densities, ionization levels, and turbulence,
conditions for planetary habitability also evolve. The heliosphere will exit
the LIC and enter a new interstellar cloud in the next 2,000 years [9], but
the precise location of the cloud boundary, and the expected effects on the
heliosphere and life within it, remain poorly constrained. The large increase
in pulsar (and other) sightlines in the next decade may sharpen our knowledge
of that local boundary.

A Galaxy of Astrospheres. The multi-phase ISM is permeated by billions of
stellar astrospheres [45], most of which host at least one exoplanet. As the
known exoplanet population continues to dramatically rise thanks to ongoing
searches with new and upcoming telescopes, including the James Webb Space
Telescope, characterization of their habitability will require consideration
of star-ISM interactions analogous to those driving the formation of the
heliosphere and the conditions for life within it. ISM conditions vary
dramatically across the Galaxy, making continued remote study of the ISM
critical to characterization of other stellar astrospheres and their
exoplanets.

##  4 The Next Decade of Opportunities

Over the next decade we will begin painting a holistic picture of the
heliosphere and its Galactic context, not only via large-scale studies of the
solar wind, but also through a global assessment of the heliosphere’s
interstellar environment and its connection to physical conditions in the
broader ISM. This work is only possible through interdisciplinary studies
enabled by a combination of spacecraft and ground-based observations, as
summarized in the table below.

## References

  * [1] Planck Collaboration et al.  “Planck 2018 results. XI. Polarized dust foregrounds”  In _A &A_ 641, 2020, pp. A11  DOI: 10.1051/0004-6361/201832618
  * [2] Douglas P. Finkbeiner  “A Full-Sky H$\alpha$ Template for Microwave Foreground Prediction”  In _ApJS_ 146.2, 2003, pp. 407–415  DOI: 10.1086/374411
  * [3] D.. Gurnett, W.. Kurth, L.. Burlaga and N.. Ness  “In Situ Observations of Interstellar Plasma with Voyager 1”  In _Science_ 341.6153, 2013, pp. 1489–1492  DOI: 10.1126/science.1241681
  * [4] L.. Burlaga, N.. Ness and E.. Stone  “Magnetic Field Observations as Voyager 1 Entered the Heliosheath Depletion Region”  In _Science_ 341.6142, 2013, pp. 147–150  DOI: 10.1126/science.1235451
  * [5] E.. Stone et al.  “Voyager 1 Observes Low-Energy Galactic Cosmic Rays in a Region Depleted of Heliospheric Ions”  In _Science_ 341.6142, 2013, pp. 150–153  DOI: 10.1126/science.1236408
  * [6] L.. Burlaga et al.  “Magnetic field and particle measurements made by Voyager 2 at and near the heliopause”  In _Nature Astronomy_ 3, 2019, pp. 1007–1012  DOI: 10.1038/s41550-019-0920-y
  * [7] Edward C. Stone, Alan C. Cummings, Bryant C. Heikkila and Nand Lal  “Cosmic ray measurements from Voyager 2 as it crossed into interstellar space”  In _Nature Astronomy_ 3, 2019, pp. 1013–1018  DOI: 10.1038/s41550-019-0928-3
  * [8] Jeffrey L. Linsky, Seth Redfield and Dennis Tilipman  “The Interface between the Outer Heliosphere and the Inner Local ISM: Morphology of the Local Interstellar Cloud, Its Hydrogen Hole, Strömgren Shells, and 60Fe Accretion”  In _ApJ_ 886.1, 2019, pp. 41  DOI: 10.3847/1538-4357/ab498a
  * [9] Jeffrey Linsky, Seth Redfield, Diana Ryder and Eberhard Moebius  “Inhomogeneity in the Local ISM and Its Relation to the Heliosphere”  In _Space Science Reviews_ 218.3, 2022, pp. 16  DOI: 10.1007/s11214-022-00884-5
  * [10] Jens Kleimann et al.  “The Structure of the Large-Scale Heliosphere as Seen by Current Models”  In _Space Science Reviews_ 218.4, 2022, pp. 36  DOI: 10.1007/s11214-022-00902-6
  * [11] D.. Gurnett and W.. Kurth  “Plasma densities near and beyond the heliopause from the Voyager 1 and 2 plasma wave instruments”  In _Nature Astronomy_ 3, 2019, pp. 1024–1028  DOI: 10.1038/s41550-019-0918-5
  * [12] John D. Richardson, John W. Belcher, Paula Garcia-Galindo and Leonard F. Burlaga  “Voyager 2 plasma observations of the heliopause and interstellar medium”  In _Nature Astronomy_ 3, 2019, pp. 1019–1023  DOI: 10.1038/s41550-019-0929-2
  * [13] D.. Gurnett et al.  “A Foreshock Model for Interstellar Shocks of Solar Origin: Voyager 1 and 2 Observations”  In _AJ_ 161.1, 2021, pp. 11  DOI: 10.3847/1538-3881/abc337
  * [14] Stella Koch Ocker et al.  “Persistent plasma waves in interstellar space detected by Voyager 1”  In _Nature Astronomy_ , 2021  DOI: 10.1038/s41550-021-01363-7
  * [15] D.. Gurnett et al.  “Origin of the Weak Plasma Emission Line Detected by Voyager 1 in the Interstellar Medium: Evidence for Suprathermal Electrons”  In _ApJ_ 921.1, 2021, pp. 62  DOI: 10.3847/1538-4357/ac1c7a
  * [16] N. Meyer-Vernet, A. Lecacheux, K. Issautier and M. Moncuquet  “Weak line discovered by Voyager 1 in the interstellar medium: Quasi-thermal noise produced by very few fast electrons”  In _A &A_ 658, 2022, pp. L12  DOI: 10.1051/0004-6361/202243030
  * [17] L.. Burlaga et al.  “Magnetic Field and Plasma Density Observations of a Pressure Front by Voyager 1 during 2020 in the Very Local Interstellar Medium”  In _ApJ_ 911.1, 2021, pp. 61  DOI: 10.3847/1538-4357/abeb6a
  * [18] Ralph L. McNutt et al.  “Interstellar probe - Destination: Universe!”  In _Acta Astronautica_ 196, 2022, pp. 13–28  DOI: 10.1016/j.actaastro.2022.04.001
  * [19] Pontus C. Brandt et al.  “Interstellar Probe: Humanity’s exploration of the Galaxy Begins”  In _Acta Astronautica_ , 2022  DOI: https://doi.org/10.1016/j.actaastro.2022.07.011
  * [20] X.. You et al.  “An Improved Solar Wind Electron Density Model for Pulsar Timing”  In _ApJ_ 671.1, 2007, pp. 907–911  DOI: 10.1086/522227
  * [21] X.. You, W.. Coles, G.. Hobbs and R.. Manchester  “Measurement of the electron density and magnetic field of the solar wind using millisecond pulsars”  In _MNRAS_ 422.2, 2012, pp. 1160–1165  DOI: 10.1111/j.1365-2966.2012.20688.x
  * [22] D.. Madison et al.  “The NANOGrav 11 yr Data Set: Solar Wind Sounding through Pulsar Timing”  In _ApJ_ 872.2, 2019, pp. 150  DOI: 10.3847/1538-4357/ab01fd
  * [23] C. Tiburzi et al.  “On the usefulness of existing solar wind models for pulsar timing corrections”  In _MNRAS_ 487.1, 2019, pp. 394–408  DOI: 10.1093/mnras/stz1278
  * [24] C. Tiburzi et al.  “The impact of solar wind variability on pulsar timing”  In _A &A_ 647, 2021, pp. A84  DOI: 10.1051/0004-6361/202039846
  * [25] Jeffrey S. Hazboun et al.  “Bayesian Solar Wind Modeling with Pulsar Timing Arrays”  In _ApJ_ 929.1, 2022, pp. 39  DOI: 10.3847/1538-4357/ac5829
  * [26] CHIME/Pulsar Collaboration et al.  “The CHIME Pulsar Project: System Overview”  In _ApJS_ 255.1, 2021, pp. 5  DOI: 10.3847/1538-4365/abfdcb
  * [27] B.. Stappers et al.  “Observing pulsars and fast transients with LOFAR”  In _A &A_ 530, 2011, pp. A80  DOI: 10.1051/0004-6361/201116681
  * [28] M.. van Haarlem et al.  “LOFAR: The LOw-Frequency ARray”  In _A &A_ 556, 2013, pp. A2  DOI: 10.1051/0004-6361/201220873
  * [29] K. Stovall et al.  “Pulsar Observations Using the First Station of the Long Wavelength Array and the LWA Pulsar Data Archive”  In _ApJ_ 808.2, 2015, pp. 156  DOI: 10.1088/0004-637X/808/2/156
  * [30] P. Kumar et al.  “Pulsar Observations at Low Frequencies: Applications to Pulsar Timing and Solar Wind Models”, 2022  DOI: 10.1093/mnras/stac316
  * [31] Laura D. Ingleby, Steven R. Spangler and Catherine A. Whiting  “Probing the Large-Scale Plasma Structure of the Solar Corona with Faraday Rotation Measurements”  In _ApJ_ 668.1, 2007, pp. 520–532  DOI: 10.1086/521140
  * [32] B. van der Holst et al.  “Alfvén Wave Solar Model (AWSoM): Coronal Heating”  In _ApJ_ 782.2, 2014, pp. 81  DOI: 10.1088/0004-637X/782/2/81
  * [33] James W. McKee, Hengrui Zhu, Daniel R. Stinebring and James M. Cordes  “Probing the Local Interstellar Medium with Scintillometry of the Bright Pulsar B1133 + 16”  In _ApJ_ 927.1, 2022, pp. 99  DOI: 10.3847/1538-4357/ac460b
  * [34] Dan R. Stinebring et al.  “A Scintillation Arc Survey of 22 Pulsars with Low to Moderate Dispersion Measures”  In _arXiv e-prints_ , 2022, pp. arXiv:2207.08756  arXiv:2207.08756 [astro-ph.HE]
  * [35] Ziwei Wu et al.  “Pulsar scintillation studies with LOFAR. I. The census”  In _A &A_ 663, 2022, pp. A116  DOI: 10.1051/0004-6361/202142980
  * [36] James M. Cordes, Barney J. Rickett, Daniel R. Stinebring and William A. Coles  “Theory of Parabolic Arcs in Interstellar Scintillation Spectra”  In _ApJ_ 637.1, 2006, pp. 346–365  DOI: 10.1086/498332
  * [37] Jeffrey L. Linsky and Seth Redfield  “Could the Local Cavity be an Irregularly Shaped Strömgren Sphere?”  In _ApJ_ 920.2, 2021, pp. 75  DOI: 10.3847/1538-4357/ac1feb
  * [38] C.. Peri et al.  “E-BOSS: an Extensive stellar BOw Shock Survey. I. Methods and first catalogue”  In _A &A_ 538, 2012, pp. A108  DOI: 10.1051/0004-6361/201118116
  * [39] C.. Peri, P. Benaglia and N.. Isequilla  “E-BOSS: An Extensive stellar BOw Shock Survey. II. Catalogue second release”  In _A &A_ 578, 2015, pp. A45  DOI: 10.1051/0004-6361/201424676
  * [40] L. Decin et al.  “The enigmatic nature of the circumstellar envelope and bow shock surrounding Betelgeuse as revealed by Herschel. I. Evidence of clumps, multiple arcs, and a linear bar-like structure”  In _A &A_ 548, 2012, pp. A113  DOI: 10.1051/0004-6361/201219792
  * [41] James M. Cordes, Roger W. Romani and Scott C. Lundgren  “The Guitar nebula: a bow shock from a slow-spin, high-velocity neutron star”  In _Nature_ 362.6416, 1993, pp. 133–135  DOI: 10.1038/362133a0
  * [42] S. Chatterjee and J.. Cordes  “Bow Shocks from Neutron Stars: Scaling Laws and Hubble Space Telescope Observations of the Guitar Nebula”  In _ApJ_ 575.1, 2002, pp. 407–418  DOI: 10.1086/341139
  * [43] Stella Koch Ocker, James M. Cordes, Shami Chatterjee and Timothy Dolch  “An In Situ Study of Turbulence near Stellar Bow Shocks”  In _ApJ_ 922.2, 2021, pp. 233  DOI: 10.3847/1538-4357/ac2b28
  * [44] Priscilla C. Frisch, Seth Redfield and Jonathan D. Slavin  “The Interstellar Medium Surrounding the Sun”  In _ARA &A_ 49.1, 2011, pp. 237–279  DOI: 10.1146/annurev-astro-081710-102613
  * [45] Brian E. Wood  “Astrospheres and Solar-like Stellar Winds”  In _Living Reviews in Solar Physics_ 1.1, 2004, pp. 2  DOI: 10.12942/lrsp-2004-2