#  Radiative Transfer For Variable 3D Atmospheres
Working Document

F. Golse, F. Hecht, O. Pironneau, D. Smets and P.-H. Tournier

###### Abstract

To study the temperature in a gas subjected to electromagnetic radiations, one
may use the Radiative Transfer equations coupled with the Navier-Stokes
equations. The problem has 7 dimensions; however with minimal simplifications
it is equivalent to a small number of integro-differential equations in 3
dimensions. We present the method and a numerical implementation using an
${\mathcal{H}}$-matrix compression scheme. The result is a very fast: 50K
physical points, all directions of radiation and 680 frequencies require less
than 5 minutes on an Apple M1 Laptop. The method is capable of handling
variable absorptioN and scattering functionS of spatial positions and
frequencies.

The implementation is done using htool111https://github.com/htool-ddm/htool, a
matrix compression library interfaced with the PDE solverfreefem++.
Applications to the temperature in the French Chamonix valley is presented at
different hours of the day with and without snow / clouds and with a variable
absorption taken from the Gemini measurements. The result is precise enough to
assert temperature differences due to increased absorption in the vibrational
frequency subrange of greenhouse gasses.

Keywords Radiative transfer, climatology, clouds, Integral equation, Navier-
Stokes equation, H-matrix, Finite Element Methods.

## Introduction

Heat transfer with radiative transfer are very important in astronomy,
combustion and climatology, to name a few. For the atmosphere the reader is
sent to [12],[9], [2], the numerically oriented book [17] and the two
mathematically oriented books [4] and [5]. Due to the considerable difference
of length scales, modelling at the level of photons is difficult to use for
the earth atmosphere in large areas. A much simpler formulation, known as the
radiative transfer equations, is based on energy conservation principles of
continuum mechanics.

These were shown to be well posed in [15] for the Radiative Transfer system
coupled with the time dependent heat equation. Existence and uniqueness was
proved in [7] and [6] for the stationary case.

On the numerical side the problem is hard because the model has 3 space
variables, two directional ones, the frequencies and time.

In 2005, K. Evans and A. Marshak wrote in chapter 4 of [12] a review of the
numerical methods available for Radiative Transfer alone. Today, judging from
[3], the situation has not changed: SHDOM (Spherical Harmonic Discrete
Ordinate Method) and Monte-Carlo are the two most popular methods. While
reviewing the current situation for the radiative transfer equations in [1] we
implemented a finite element version of SHDOM and found that the method was
incapable, unless a huge number of degree of freedom is used, of giving
results with the accuracy needed to differentiate between small variations of
the absorption coefficient.

One may resort to approximations. In nuclear physics and astronomy, for the
numerical simulations in the 1960 a constant absorption coefficient was used,
the so called Grey Model. In climatology, the grey model cannot explain the
role of the greenhouse gasses.

The stratified approximation assumes that the surface which receives the
radiations is flat and the source is far. Thus only one space variable is
retained. For such cases, an integral formulation, probably due to
Chandrasekhar [4], turns out to be much more precise and also computationally
much cheaper. A fixed-point iteration scheme on this nonlinear integral
formulation, known in the radiative transfer community as “iterations on the
sources” was shown to be monotone in [14], a property which seems to have
escaped earlier studies. Finally in [7],and [8] the method was extended to
include the temperature equation of the fluid and also to handle Rayleigh
scattering while retaining monotonicity. In [6] F. Golse observed that the
integral formulation exists also in 3D. A new existence result was proved by
using this formulation and it was numerically tested. However the computing
cost was $O(N^{2})$, $N$ being the number of vertices of the triangulation of
the physical domain.

The purpose of this article is to present an $O(N\ln N)$ implementation which
uses compressed ${\mathcal{H}}$-matrices and views the computation of
integrals as a matrix-vector product in the Finite Element discrete space.
Thence P.H. Tournier, who is a co-author of the library htool for integral
equations with ${\mathcal{H}}$-matrices, wrote the necessary modifications for
Radiative Transfer, namely the handling of mixed matrices with one vertex on
the boundary and the other one inside the domain.

The programming is a joint work of F. Hecht, O. Pironneau and P.H. Tournier,
done using the high level PDE solver freefem++ [10].

The method is tested on a $35\times 35$ km area and the first 10km of the
atmosphere above it. The center of the domain is the city of Chamonix in the
French Alps, provided by D. Smets who has constructed high precision maps and
an automatic triangular mesh generator for many places in the world for the
Microsoft flight simulator. Hence this Radiative-Heat-Transfer code can be
used easily for any other terrain. Note that such a study would be a challenge
for Cartesian meshes because some of the mountains like the Mont-Blanc are
above 4000m and some slopes are almost vertical.

The absorption coefficient is a function of space and frequencies. It could
also be a function of temperature if another fixed point loop is added. No
regularity is assumed beside positivity and boundedness. Similarly, the
scattering variable coefficient need only to be in (0,1).

Real life data are used; the scattering part is not yet implemented As shown
in [FGOP2], it is not so difficult to do, but it adds complexity to an already
complex topic. Numerically, the methods is very fast: it can handle a few
hundred thousand physical points plus all directions plus 680 frequency points
on a PC (and a few millions on a supercomputer) within minutes.

Note that the authors have no competence for climate modelling, so the results
are only briefly commented: the purpose of this study is to validate the
numerical tool.

##  1 The Mathematical Problem

To find the temperature $T$ in an incompressible fluid exposed to
electromagnetic radiations, it is necessary to solve the Navier-Stokes
equations coupled with the Radiative Transfer equations. It is a system of
partial differential equations formulated in terms of the fluid velocity ${\bf
u}$, its pressure $p$, its density $\rho$, and temperature $T$; all are
functions of time $t$ and position ${\bf x}$ in the physical domain $\Omega$.
It involves also the light intensity field ${I_{\nu}}({\bf
x},t,{\bm{\omega}},t)$ for each frequency $\nu$ in each direction
${\bm{\omega}}$:

Given ${I_{\nu}},T,{\bf u},\rho$ at $t=0$, find ${I_{\nu}},T,{\bf u},\rho$,
s.t. $\\{{\bf
x},{\bm{\omega}},t,\nu\\}\in\Omega\times{\mathbb{S}}_{2}\times(0,\bar{T})\times{\mathbb{R}}^{+}$,

$\displaystyle\frac{1}{c}\partial_{t}{I_{\nu}}+{\bm{\omega}}\cdot\bm{\nabla}{I_{\nu}}+\rho\tilde{\kappa}a\left[{I_{\nu}}-{\frac{1}{4\pi}\int_{{\mathbb{S}}^{2}}}p_{\nu}({\bm{\omega}},{\bm{\omega}}^{\prime}){I_{\nu}}({\bm{\omega}}^{\prime}){d}\omega^{\prime}\right]$
(1) $\displaystyle\hskip
170.71652pt=\rho\tilde{\kappa}(1-a)[B_{\nu}(T)-{I_{\nu}}],$
$\displaystyle\rho(\partial_{t}T+{\bf u}\cdot\nabla
T)-\nabla\cdot(\rho\kappa_{T}\nabla
T)+A\bm{\nabla}\cdot\int_{0}^{\infty}\frac{1}{4\pi}\int_{{\mathbb{S}}^{2}}{I_{\nu}}({\bm{\omega}}){\bm{\omega}}\mathrm{d}\omega\mathrm{d}\nu=0$
$\displaystyle\partial_{t}{\bf u}+{\bf u}\cdot\nabla{\bf
u}-\frac{\mu_{F}}{\rho}\Delta{\bf u}+\frac{1}{\rho}\nabla p={\bf
g}(T),\quad\nabla\cdot{\bf u}=0,\quad\partial_{t}\rho+\nabla\cdot(\rho{\bf
u})=0,$

where $A=4\pi/c_{P}$, $c_{P}$ is the specific heat of the medium at constant
pressure, $\SS$ is the unit sphere, $\nabla,\Delta$ are with respect to ${\bf
x}$, $\displaystyle B_{\nu}(T)=\frac{2\hbar\nu^{3}}{c^{2}[{\rm
e}^{\frac{\hbar\nu}{kT}}-1]}$, is the Planck function, $\hbar,c,k$ are the
Planck constant, the speed of light in the medium and the Boltzmann constant.
The above holds only if $c_{P}$ is constant in $\Omega$, which is the case in
moderate size domains.

###  1.1 Adimensionalization

Recall that $c=2.99\cdot 10^{9}m\cdot s^{-1},~{}\hbar=6.63\cdot
10^{-34}J.s,~{}k=1.38\cdot 10^{-23}J\cdot K^{-1},$

$\kappa_{T}=2\cdot 10^{-5}m^{2}\cdot s^{-1},\quad
c_{P}=1.1~{}J.s\cdot(g.K)^{-1},\quad\rho\sim\bar{\rho}:=1.3g\cdot
m^{-3},\hbox{ in air at $20^{o}$C}.$

The frequencies of interest are in the order of $\nu_{0}=10^{14}s^{-1}$, so we
define $B_{0}=\frac{2\hbar\nu_{0}^{3}}{c^{2}}$, $T_{0}=\frac{\hbar\nu_{0}}{k}$
and work with $I_{\nu}/B_{0}$, $T/T_{0}$ and $\nu/\nu_{0}$. Then (LABEL:onea)
holds with these new unknowns with $\displaystyle B_{\nu}(T)={\nu^{3}}/({\rm
e}^{\frac{\nu}{T}}-1)$, $c_{P}=1$ and

$B_{0}=1.47~{}J\cdot m^{-2}\quad
T_{0}=4789K,\quad\frac{A}{\bar{\rho}}=\frac{4\pi
B_{0}\nu_{0}}{c_{P}\bar{\rho}T_{0}}=2.70\cdot 10^{11}~{}m\cdot s^{-1}.$ (2)

###### Notation 1

If $\bar{\rho}$ is a reference density for the gas,
$\kappa:=\bar{\rho}\tilde{\kappa}$ is the percentage of light absorbed per
unit length. It has the dimension $m^{-1}$, thence, the notation

$\kappa({\bf x},\nu):=\rho({\bf x})\tilde{\kappa}_{\nu}({\bf x}).$

Later we may assume at times that $\tilde{\kappa}_{\nu}({\bf x})$ is
independent of ${\bf x}$ or is the sum of products of ${\bf x}$-functions with
$\nu$-functions.

The scattering albedo is $a\in(0,1)$ and
$\frac{1}{4\pi}p_{\nu}({\bm{\omega}},{\bm{\omega}}^{\prime})$ is the
“probability” that a ray in direction ${\bm{\omega}}^{\prime}$ scatters in
direction ${\bm{\omega}}$. Both $\kappa$ and $a$ are functions of $\nu$. We
further assume that $\nu\mapsto\kappa$ and $\nu\mapsto a$ are continuous
positive functions, satisfying, for some positive constants $a_{M}$ and
$\kappa_{m}<\kappa_{M}$:

$0\leq\kappa_{m}\leq\kappa\leq\kappa_{M}\,,\qquad 0\leq a\leq
a_{M}<1\,,\qquad\forall\nu>0\,.$

The viscosity of air $\mu_{F}$ is $18\cdot 10^{-6}$Pa.s; ${\bf g}(T)$, the
Boussinesq term, is a vector valued function of the temperature $T$.

We require $\Omega$ to be an open bounded subset of $\mathbf{R}^{3}$ with
$C^{1}$ boundary. We denote by ${\bf n}({\bf x})$ the outward unit normal to
$\Gamma:={\partial}\Omega$ at ${\bf x}$.

Boundary conditions are needed, for example, initial values for
$I_{\nu},T,{\bf u},\rho$, ${\bf u}$ given on $\Gamma$ and $\rho$ given at
${\bf x}\in\Gamma$ where ${\bf u}\cdot{\bf n}<0$ and

$I_{\nu}({\bf x},{\bm{\omega}})\\!=\\!Q_{\nu}({\bf
x},{\bm{\omega}})\,,\quad{\bm{\omega}}\cdot{\bf
n}<0\,,\qquad\frac{{\partial}T}{{\partial}n}\Big{|}_{\Gamma}=0\,.$

####  1.1.1 Sunlight

When $Q_{\nu}({\bf x},{\bm{\omega}})$ is due to sunlight, of power $Q^{0}$
coming from the direction ${\bm{\omega}}_{s}$, then the surface which receives
the light re-emits it with (Lambertian reflexion),

$Q_{\nu}({\bf x},{\bm{\omega}})=Q_{\nu}({\bf x})\cdot[{\bm{\omega}}\cdot{\bf
n}({\bf x})]_{-}\hbox{ with }Q_{\nu}({\bf
x}):=Q^{0}[{\bm{\omega}}_{s}\cdot{\bf n}({\bf x})]_{+}B_{\nu}(T_{s}),$

the Planck function with the Sun’s temperature $T_{s}$.

###  1.2 Simplifications

The angular average radiative intensity plays an important role in this
article:

$J_{\nu}({\bf
x},t)=\frac{1}{4\pi}\int_{{\mathbb{S}}^{2}}{I_{\nu}}({\bm{\omega}}){d}\omega.$

Averaging the first equation in (LABEL:onea) and neglecting the time dependent
term because $c>>1$, gives :

$\bm{\nabla}\cdot\int_{0}^{\infty}{\int_{{\mathbb{S}}^{2}}}{I_{\nu}}({\bm{\omega}}){\bm{\omega}}{d}\omega{d}\nu=4\pi\int_{0}^{\infty}\rho\tilde{\kappa}(1-a)\left(B_{\nu}(T)-J_{\nu}\right){d}\nu\,,$
(3)

which, in turn, leads to a simpler form ( see (4)) for the second equation in
(LABEL:onea).

####  1.2.1 Isotropic Scattering

Assume isotropic scattering:
$p_{\nu}({\bm{\omega}},{\bm{\omega}}^{\prime})\equiv 1$. Let us neglect the
variations of $\rho$ in the diffusion term of the temperature equation. Let us
neglect $\frac{1}{c}\partial_{t}I$ because $c$ is very large and let us assume
that the Boussinesq term is small. Then the Navier-Stokes equations are
decoupled and can be solved before hand. Then in the Radiative transfer and
temperature equations $\rho$, ${\bf u}$ are given. So (LABEL:onea) becomes the
following system for $I_{\nu},J_{\nu},T$:

$\left\\{\begin{aligned} {}&{\bm{\omega}}\cdot\nabla
I_{\nu}+{\rho\tilde{\kappa}}I_{\nu}={\rho\tilde{\kappa}}(1-a)B_{\nu}(T)+{\rho\tilde{\kappa}}aJ_{\nu}\,,\quad
J_{\nu}:=\tfrac{1}{4\pi}\int_{\SS}I_{\nu}\mathrm{d}{\bm{\omega}}\,,\\\
&\partial_{t}T+{\bf u}\cdot\nabla T-\nabla\cdot(\kappa_{T}\nabla
T)=A\int_{0}^{\infty}{\tilde{\kappa}}(1-a)(J_{\nu}-B_{\nu}(T))\mathrm{d}\nu\,,\\\
&I_{\nu}({\bf x},{\bm{\omega}})\\!=\\!Q_{\nu}({\bf
x},{\bm{\omega}})\,,\quad{\bm{\omega}}\cdot{\bf n}<0\,,\,\,{\bf
x}\in{\partial}\Omega\,,\qquad\frac{{\partial}T}{{\partial}n}\Big{|}_{{\partial}\Omega}=0\,,\quad
T\big{|}_{t=0}=T_{0}\,.\end{aligned}\right.$ (4)

####  1.2.2 Small Thermal Diffusion and Convection

When the convective velocity and $\kappa_{T}$ are small compared to $A$, the
temperature equation reaches a stationary states and simplifies to an integral
relation between $J$ and $B_{\nu}(T)$ and the system becomes:

$\displaystyle{\bm{\omega}}\cdot\nabla
I_{\nu}+{\rho\tilde{\kappa}}I_{\nu}={\rho\tilde{\kappa}}(1-a)B_{\nu}(T)+{\rho\tilde{\kappa}}aJ_{\nu}\,,\quad
J_{\nu}:=\tfrac{1}{4\pi}\int_{\SS}I_{\nu}\mathrm{d}{\bm{\omega}}\,,$ (5)
$\displaystyle\int_{0}^{\infty}{\rho\tilde{\kappa}}(1-a)(J_{\nu}-B_{\nu}(T))\mathrm{d}\nu=0\,,$
(6) $\displaystyle I_{\nu}({\bf x},{\bm{\omega}})=Q_{\nu}({\bf
x},{\bm{\omega}})\,\hbox{ on }\Sigma:=\\{({\bf
x},{\bm{\omega}})\in\Gamma\times\SS~{}:~{}{\bm{\omega}}({\bf x})\cdot{\bf
n}({\bf x})<0\,\\}\,.$ (7)

The system is 6 dimensional in $\left[{\bf
x}=[x,y,z],{\bm{\omega}}=[\omega_{1},\omega_{2},\omega_{3}~{}:~{}|{\bm{\omega}}|=1],\nu\right]$.

####  1.2.3 The stratified case

When $\Omega$ is thin in a direction $Ox$ and the physical data depend mostly
on $x$, the equations becomes almost one dimensional in $x$. The domain
$\Omega$ is locally $(0,H)\times{\mathbb{R}}^{2}$. Assuming that sunlight at
$x=+\infty$ crosses the atmosphere unaffected the equations reduce to

$\displaystyle\mu\partial_{x}I_{\nu}+{\rho\tilde{\kappa}}I_{\nu}={\rho\tilde{\kappa}}(1-a)B_{\nu}(T)+{\rho\tilde{\kappa}}aJ_{\nu}\,,\quad
J_{\nu}:=\tfrac{1}{2}\int_{0}^{1}I_{\nu}\mathrm{d}\mu$ (8)
$\displaystyle\int_{0}^{\infty}{\rho\tilde{\kappa}}(1-a)(J_{\nu}-B_{\nu}(T))\mathrm{d}\nu=0\,,$
$\displaystyle I_{\nu}(0,\mu)=Q_{\nu}({\bf x},\mu),~{}~{} \mu<0,\quad
I_{\nu}(H,\mu)=0,~{}~{}\mu>0.$

with $\Gamma$ the plane $x=0$, $Q_{\nu}=Q^{0}B_{\nu}(T_{s})\mu$,
$\mu:=\cos\alpha$ and $\alpha$ the angle of the ray ${\bm{\omega}}$ with $Ox$.

##  2 The Method of Characteristics

Let us derive a closed form solution of (5) by computing the long time
solution of

$\partial_{\tau}I({\bf x},\tau)+{\bm{\omega}}\nabla_{\bf x}I({\bf
x},\tau)+\kappa({\bf x})I({\bf x},\tau)=S({\bf x}).$ (9)

While it plays no part, notice that $\tau=t/c$. Let $K({\bf
x},{\bm{\omega}},\tau)=\int_{0}^{\tau}\kappa({\bf
x}+{\bm{\omega}}s)\mathrm{d}s$. Then

$\frac{\mathrm{d}}{\mathrm{d}\tau}\left(I({\bf x}+{\bm{\omega}}\tau,\tau){\rm
e}^{K({\bf x},{\bm{\omega}},\tau)}\right)={\rm e}^{K({\bf
x},{\bm{\omega}},\tau)}\left[\partial_{\tau}I+{\bm{\omega}}\nabla_{\bf
x}I+\kappa({\bf x})I\right]_{|_{{\bf
x}+{\bm{\omega}}\tau,{\bm{\omega}},\tau}}$

Consequently

$I({\bf x}+{\bm{\omega}}\tau,\tau){\rm e}^{K({\bf
x},{\bm{\omega}},\tau)}=I({\bf x},0)+\int_{0}^{\tau}{\rm e}^{K({\bf
x},{\bm{\omega}},s)}S({\bf x}+{\bm{\omega}}s)\mathrm{d}s.$

Denote ${\bf x}^{\prime}={\bf x}+{\bm{\omega}}\tau$; the above is also

$I({\bf x}^{\prime},\tau){\rm e}^{K({\bf
x}^{\prime}-{\bm{\omega}}\tau,{\bm{\omega}},\tau)}=I({\bf
x}^{\prime}-{\bm{\omega}}\tau,0)+\int_{0}^{\tau}{\rm e}^{K({\bf
x}^{\prime}-{\bm{\omega}}\tau,{\bm{\omega}},s)}S({\bf
x}+{\bm{\omega}}s)\mathrm{d}s.$

Denote by $\tau_{{\bf x}^{\prime},{\bm{\omega}}}$ the value of $\tau$ which
brings ${\bf x}^{\prime}-{\bm{\omega}}\tau$ to the boundary:

${\bf x}^{\prime}-{\bm{\omega}}\tau_{{\bf
x}^{\prime},{\bm{\omega}}}\in\Sigma,\hbox{ and denote }{\bf x}_{\Sigma}({\bf
x}^{\prime},{\bm{\omega}}):={\bf x}^{\prime}-{\bm{\omega}}\tau_{{\bf
x}^{\prime},{\bm{\omega}}}.$

Then the long time solution of (9) is

$\displaystyle I({\bf x}^{\prime})$ $\displaystyle=I({\bf x}_{\Sigma}({\bf
x}^{\prime},{\bm{\omega}})){\rm e}^{-K({\bf x}_{\Sigma}({\bf
x}^{\prime},{\bm{\omega}}),{\bm{\omega}},\tau)}$ $\displaystyle\hskip
56.9055pt+\int_{0}^{\tau_{{\bf x}^{\prime},{\bm{\omega}}}}{\rm e}^{K({\bf
x}^{\prime}-{\bm{\omega}}s,{\bm{\omega}},s)-K({\bf x}_{\Sigma}({\bf
x}^{\prime},{\bm{\omega}}),{\bm{\omega}},\tau_{{\bf
x}^{\prime},{\bm{\omega}}})}S({\bf x}^{\prime}-{\bm{\omega}}(\tau_{{\bf
x}^{\prime},{\bm{\omega}}}-s))\mathrm{d}s.$

It can also be written as:

$\displaystyle I({\bf x},{\bm{\omega}})$ $\displaystyle=I({\bf
x}_{\Sigma}({\bf x},{\bm{\omega}})){\rm e}^{-\int_{0}^{\tau_{{\bf
x},{\bm{\omega}}}}\kappa({\bf
x}-{\bm{\omega}}s)\mathrm{d}s}+\int_{0}^{\tau_{{\bf x},{\bm{\omega}}}}{\rm
e}^{-\int_{0}^{s}\kappa({\bf
x}-{\bm{\omega}}s^{\prime})\mathrm{d}s^{\prime}}S({\bf
x}-{\bm{\omega}}s)\mathrm{d}s$ (10) $\displaystyle=I({\bf x}_{\Sigma}({\bf
x},{\bm{\omega}})){\rm e}^{-\int_{[{\bf x},{\bf
x}_{\Sigma}]}\kappa}+\int_{[{\bf x},{\bf x}_{\Sigma}]}{\rm e}^{-\int_{[{\bf
x},{\bf x}^{\prime}]}\kappa}S({\bf x}^{\prime})\mathrm{d}x^{\prime}\,.$

Averaging in ${\bm{\omega}}$ and using (7), and knowing that the solid angle
integral of $f(x-{\bm{\omega}}s)\mathrm{d}\omega$ is the surface integral on
$\Gamma$ of $f({\bf x}^{\prime}-{\bf x})\cdot{\bf n}({\bf x}^{\prime})|{\bf
x}^{\prime}-{\bf x}|^{-3}\mathrm{d}\Gamma(x^{\prime})$ leads to :

$\displaystyle J_{\nu}({\bf x})$
$\displaystyle:=\frac{1}{4\pi}\int_{\SS}I({\bf
x},{\bm{\omega}})\mathrm{d}\omega=\frac{1}{4\pi}\int_{\SS}I({\bf
x}_{\Sigma}({\bf x},{\bm{\omega}})){\rm e}^{-\int_{0}^{\tau_{{\bf
x},{\bm{\omega}}}}\kappa({\bf x}-{\bm{\omega}}s)\mathrm{d}s}\mathrm{d}\omega$
(11) $\displaystyle+\frac{1}{4\pi}\int_{\SS}\int_{0}^{\tau_{{\bf
x},{\bm{\omega}}}}{\rm e}^{-\int_{0}^{s}\kappa({\bf
x}-{\bm{\omega}}s^{\prime})\mathrm{d}s^{\prime}}S({\bf
x}-{\bm{\omega}}s)\mathrm{d}s\mathrm{d}\omega$ $\displaystyle=S^{E}_{\nu}({\bf
x})+{\mathcal{J}}[S]({\bf x}).$

with the notations, $z_{-}:=-\min(z,0)$ and

$\displaystyle S^{E}_{\nu}({\bf x}):==\frac{1}{4\pi}\int_{\Gamma}Q_{\nu}({\bf
x},\frac{{\bf x}^{\prime}-{\bf x}}{|{\bf x}^{\prime}-{\bf x}|})\frac{[({\bf
x}^{\prime}-{\bf x})\cdot{\bf n}({\bf x}^{\prime})]_{-}}{|{\bf
x}^{\prime}-{\bf x}|^{3}}{\rm e}^{-\int_{[{\bf x},{\bf
x}^{\prime}]}\kappa}\mathrm{d}\Gamma({\bf x}^{\prime})$ (12)
$\displaystyle{\mathcal{J}}[S]({\bf x}):=\frac{1}{4\pi}\int_{\Omega}S({\bf
x}^{\prime})\frac{{\rm e}^{-\int_{[{\bf x},{\bf x}^{\prime}]}\kappa}}{|{\bf
x}^{\prime}-{\bf x}|^{2}}\mathrm{d}x^{\prime}.$

where ${\bf n}({\bf x}^{\prime})$ is the outer normal at ${\bf
x}^{\prime}\in\partial\Omega$. The last line is true only if $\Omega$ is
convex but it can be used in the general case with the following definition of
$\tilde{\Omega}$ and $\bar{\kappa}$.

${\bf x}^{i}$${\bf x}^{\prime\prime}$${\bf x}^{\prime}$ Figure 1: Consider
$\Omega$, the domain between the 2 spheres. The segment $[{\bf x}^{i},{\bf
x}^{\prime}]$ participates to the radiative intensity in ${\bf x}^{i}$.
However when the segment $[{\bf x}^{i},{\bf x}^{\prime\prime}]$ crosses the
inner sphere, the light scattered from ${\bf x}^{\prime\prime}$ does not
contribute to the light intensity in ${\bf x}^{i}$.

###### Proposition 1

Let $\bar{\Omega}\supset\Omega$ be a convex extension of $\Omega$. Let
$\bar{\kappa}$ be an extension of $\kappa$ in $\bar{\Omega}$ such that

$\bar{\kappa}({\bf x})=\kappa(vx),~{}~{}{\bf
x}\in\Omega,\quad\bar{\kappa}({\bf x})=+\infty,~{}~{}{\bf
x}\in\bar{\Omega}\backslash\Omega.$

Then

$\int_{\SS}\int_{0}^{\tau_{{\bf x},{\bm{\omega}}}}{\rm
e}^{-\int_{0}^{s}\kappa({\bf
x}-{\bm{\omega}}s^{\prime})\mathrm{d}s^{\prime}}S({\bf
x}-{\bm{\omega}}s)\mathrm{d}s\mathrm{d}\omega=\int_{\bar{\Omega}}S({\bf
x}^{\prime})\frac{{\rm e}^{-\int_{[{\bf x},{\bf
x}^{\prime}]}\bar{\kappa}}}{|{\bf x}-{\bf
x}^{\prime}|^{2}}\mathrm{d}x^{\prime}.$

Proof : If $\Omega$ is convex then $\\{{\bf
x}-{\bm{\omega}}s~{}:~{}s\in(0,\tau_{{\bf
x},{\bm{\omega}}}),{\bm{\omega}}\in\SS\\}=\Omega$. Refering to figure 1, if
$\Omega$ is not convex, then the part in $\bar{\Omega}\backslash\Omega$ of the
integral on the right is zero because ${\bf
x}^{\prime}\in\bar{\Omega}\backslash\Omega$ implies that a portion of the
integral in the exponential will be on $[{\bf x}^{\prime\prime},{\bf
x}^{\prime}]$, for some ${\bf
x}^{\prime\prime}\in\bar{\Omega}\backslash\Omega$. This portion will
contribute to the integral with $\bar{\kappa}=+\infty$ and ${\rm
e}^{-\infty}=0$. $\square$

####  2.0.1 The Stratified Grey Case

The following will be helpful to assert the precision of the general algorithm
2.1 below.

Recall the Stefan-Blotzmann relation: $\int_{0}^{\infty}B_{\nu}(T)=\sigma
T^{4}$ with $\sigma=\pi^{4}/15$, the Stefan constant. If
$\tilde{\kappa}_{\nu}$ is constant (grey case), then (8) integrated in $\nu$
yields an equation for $\bar{I}=\int_{0}^{\infty}I_{\nu}\mathrm{d}\nu$:

$\displaystyle\mu\partial_{x}\bar{I}+\rho\tilde{\kappa}\bar{I}=\rho\tilde{\kappa}\sigma
T^{4}(x)\,,\quad\sigma
T^{4}(x)=\frac{1}{2}\int_{-1}^{1}\bar{I}(x,\mu)\mathrm{d}\mu,\quad x\in(0,H)$
(13)

with boundary conditions $\bar{I}(0,\mu)|_{\mu>0}=\mu Q^{0}\sigma T_{s}^{4}$
and $\bar{I}(H,\mu)|_{\mu<0}=0$.

Even though one can find directly a solution in integral form, when
$\rho\tilde{\kappa}$ is constant the above framework yields

$\bar{S}^{E}({\bf x})=\frac{Q^{0}\sigma
T_{s}^{4}}{4\pi}\int_{{\mathbb{R}}^{2}}\left(\frac{\cos^{2}\theta}{x}\right)^{2}{\rm
e}^{-\frac{\kappa x}{\cos\theta}}\mathrm{d}\Gamma(y,z)=\frac{Q^{0}\sigma
T_{s}^{4}}{2}{\bf E}_{3}(\kappa x)$ (14)

where ${\bf E}_{3}$ is the third exponential integral. Similarly (see [8])

${\cal J}[S](x)=\kappa\int_{0}^{H}S(x^{\prime}){\bf
E}_{1}(\kappa[x^{\prime}-x|)\mathrm{d}x^{\prime}.$ (15)

This approximation can be extended to the case $\kappa$ function of $x$ by a
change of coordinate $x\to\tau$ with
$\kappa\frac{{\partial}}{{\partial}\tau}=\frac{{\partial}}{{\partial}x}$.

###  2.1 Algorithm

System (4) can be solved by the following iterative scheme.

  1. 1.

Start from $T^{0}\equiv 0$ and $J_{\nu}^{0}=S^{E}_{\nu}$, given by (12)

  2. 2.

FOR $n=0,1,\dots,N-1$

    1. (a)

Knowing $T^{n}({\bf x},t)$ and ${J_{\nu}}^{n}({\bf x},t)$, compute by (11)

${J_{\nu}}^{n+1}=S_{\nu}^{E}+\mathcal{J}[aJ_{\nu}^{n}+(1-a)B_{\nu}(T^{n})],\;\nu\in(0,\infty),~{}{\bf
x}\in\Omega,~{}t<T.$

    2. (b)

Define $T^{n+1}$ the solution of the semilinear drift-diffusion equation

$\left\\{\begin{aligned} &{}\rho c_{V}({\partial}_{t}T^{n+1}+{\bf
u}\cdot\nabla_{x}T^{n+1})-\nabla_{x}\cdot(\rho
c_{P}\kappa_{T}\nabla_{x}T^{n+1})\\\ &\hskip
56.9055pt+\mathcal{B}(T^{n+1})=4\pi\int_{0}^{\infty}\kappa(1-a)J_{\nu}^{n+1}\mathrm{d}\nu\,,\\\
&T^{n+1}\big{|}_{t=0}=T_{0}\,,\qquad\frac{{\partial}T^{n+1}}{{\partial}n}\Big{|}_{{\partial}\Omega}=0\,,\quad{\bf
x}\in\Omega\,,\,\,t>0\,,\end{aligned}\right.$

where $\rho$ and ${\bf u}$ are given by the Navier-Stokes equations and

$\displaystyle\mathcal{B}(T):=4\pi\int_{0}^{\infty}\kappa(1-a)B_{\nu}(\min(T_{+},T_{M}))\mathrm{d}\nu\,.$

###### Remark 1

When ${\bf u}$ and $\kappa_{T}$ are negligeable, Step (b) becomes: find
$T^{n+1}$ such that

$\int_{0}^{\infty}\kappa(1-a)B_{\nu}(T_{+}^{n+1})\mathrm{d}\nu=\int_{0}^{\infty}\kappa(1-a)J_{\nu}^{n+1}\mathrm{d}\nu\,.$
(16)

This equation has a unique solution because $T\mapsto B_{\nu}(T)$ is
increasing. As for existence, observe that the left hand-side is continuous in
$T_{+}$, vanishes for $T_{+}=0$, and tends to $+\infty$ when
$T_{+}\to+\infty$. For a computer solution we may use a Newton method because
the Hessian is positive and never vanishes.

##  3 Implementation with ${\mathcal{H}}$-Matrices

Note that the computations of $S_{\nu}^{E}$ and ${\cal J}[S]$ involve
convolution like operators. Indeed, with $\Omega$ convex bounded and $a=0$,
$S^{E}$ is given by (12) and

$\displaystyle{\mathcal{J}}[S_{\nu}]({\bf x})$
$\displaystyle=\tfrac{1}{4\pi}\int_{\Omega}S_{\nu}({\bf
x}^{\prime})\kappa({\bf x}^{\prime})\frac{{\rm e}^{-\int_{[{\bf x},{\bf
x}^{\prime}]}\kappa}}{|{\bf x}-{\bf x}^{\prime}|^{2}}\mathrm{d}x^{\prime}\,.$
(17)

The domain is discretized by the vertices $\\{{\bf x}^{i}\\}_{1}^{N}$ of a
tetraedral mesh. The integrals are computed with a quadrature rule using
quadrature points inside the tetraedras, typically 25 points when ${\bf
x}^{i}$ is near to the tetraedra of the integral and 5 points otherwise ; $J$
is approximated by its $P^{1}$ interpolation on the mesh:

$J({\bf x})=\sum_{1}^{N}J_{j}{\hat{w}}^{j}({\bf x})~{}\text{ where
${\hat{w}}^{j}$ is the $P^{1}$- Finite Element hat function of vertex ${\bf
x}^{j}$}.$

Then Step (a) of Algorithm 2.1 becomes:

$\displaystyle S_{\nu,j}=aJ_{\nu,j}+(1-a)B_{\nu}(T_{j})$ (18) $\displaystyle
J_{\nu,i}=S^{E}_{\nu,i}+\sum_{j}G_{\kappa}^{ij}S_{\nu,j}\quad\text{ where
}G_{\kappa}^{ij}=\frac{1}{4\pi}\int_{\Omega}\kappa\frac{{\rm e}^{-\int_{[{\bf
x}^{i}-{\bf x}^{\prime}]}\kappa}}{|{\bf x}^{i}-{\bf
x}^{\prime}|^{2}}{\hat{w}}^{j}({\bf x}^{\prime})\mathrm{d}x^{\prime}$

The matrix ${\bf G}$ can be compressed with the ${\mathcal{H}}$-matrix method
so that the multiplication ${\bf G\cdot S}$ has complexity $O(N\ln N)$ for
each $\nu$. The method works best when the kernel of $G$ in the integral
decays exponentially with the distance between ${\bf x}^{i}$ and ${\bf
x}^{\prime}$. The ${\mathcal{H}}$-matrix approximation views $\\{G^{ij}\\}$ as
a hierarchical tree of square blocks; The blocks correspond to interaction
between clusters of points near ${\bf x}^{j}$ and near ${\bf x}^{\prime}$.
Far-field interaction blocks can be approximated by a low rank matrix because
their singular value decompositions (SVD) have fast decaying singular values.
We use partially pivoted adaptive cross-approximation [16] to approach the
first terms of the SVD of the blocks, because only r-rows times r-columns
columns are needed instead of the whole block, where r is the rank of the
approximation. The rank is a function of a user defined parameter $\epsilon$
in connection with the relative Frobenius norm error. Another criteria must be
met: if $R_{1}$ (resp $R_{2}$) is the radius of a cluster of points centered
at ${\bf x}_{1}$ (resp ${\bf x}_{2}$), then one goes down the hierarchical
tree until the corresponding block satisfies $\max(R_{1},R_{2})<\eta|{\bf
x}_{1}-{\bf x}_{2}|$ where $\eta$ is a user defined parameter. If the en of
the tree is reached, the bloc is not compacted and it is displayed in red on
figure 10.

Then, to compute the integral $\int_{0}^{\infty}\kappa J_{\nu,i}\mathrm{d}\nu$
a rectangular quadrature at points $\nu_{k}$ is used.

For example when $a=0$,

$\int_{0}^{\infty}\kappa
J_{\nu,i}\mathrm{d}\nu\approx\sum_{k}S^{E}_{\nu_{k},i}(\nu_{k}-\nu_{k-1})+\sum_{k}\sum_{j}\kappa_{\nu_{k}}G_{\kappa_{\nu_{k}}}^{ij}B_{\nu_{k}}({T_{j}})(\nu_{k}-\nu_{k-1})$
(19)

The same decomposition into product of $\cal H$-matrices with vectors
$\\{Q_{\nu_{k}}({\bf x}^{i})\\}_{i=1}^{N}$ can be applied to compute
$S^{E}_{\nu_{k}}$:

$S^{E}_{\nu_{k},i}=\sum_{1}^{N}S^{ij}_{\kappa_{\nu_{k}}}Q_{\nu_{k}}({\bf
x}^{j}),\quad S^{ij}_{\kappa}=\int_{\Gamma}\hat{w}^{j}{\rm e}^{-\int_{[{\bf
x}^{i},{\bf x}^{\prime}]}\kappa}\left(\frac{[({\bf x}^{\prime}-{\bf
x}^{i})\cdot{\bf n}(x^{\prime})]^{+}}{|{\bf x}^{\prime}-{\bf
x}^{i}|^{2}}\right)^{2}\mathrm{d}s(x^{\prime})$ (20)

###  3.1 Lebesgue Integrals

We will use the Gemini measurements to define $\kappa_{\nu}$ (figure 2). To
represent such a function we require 683 $\nu$-points. As it is, the numerical
method requires an ${\mathcal{H}}$-matrix for each $\nu$, but it is not
feasible to define $2\times 683$ ${\mathcal{H}}$-matrices!

$0$$5$$10$$15$$20$$25$$30$$0$$0.2$$0.4$$0.6$$0.8$$1$$1.2$wave lengthabsorption
coefficientapproximated $\kappa_{\nu}$$\kappa_{\nu}$

Figure 2:  Absorption $\kappa_{\nu}$ versus wavelength ($c/\nu$) in the range
$(0.1,26)$ read from the Gemini data site[11]. In the ${\mathcal{H}}$\- step
of the algorithm it is approximated by $0.01+\texttt{round}(10\kappa)/10$,
shown here as “approximated $\kappa_{\nu}$”. Right: the x-dependency of
$\kappa$, the topography and the mesh. The Mont-Blanc is in the bottom left
part and Chamonix is in the center. Colors are levels of the ${\bf x}$
dependence of $\kappa$, $0.5(1-x/2)$; the red part is the intersection of the
cloud with the mountains; the cloud is a cylinder circumscribing the red
parts, centered at Chamonix and between altitude 2000m and 8000m; it is used
only in one simulation.

Observe that the ${\mathcal{H}}$-matrices ${\mathcal{H},H}^{E}$, depend on
$\nu$ directly but through $\kappa_{\nu}$. It is an opportunity to reduce the
number of ${\mathcal{H}}$-matrices down to the number of different values of
$\kappa_{\nu}$. The the integrals with respect to $\nu$ are evaluated much
like a Lebesgue integral.

Suppose $\kappa$ takes only 2 values $\kappa_{1}$ and $\kappa_{2}$, then,
according to (20), only 2 functions $S_{\nu}^{E}$ are need,
$S^{E}_{1},S^{E}_{2}$, and similarly for $G_{\kappa}$. Now (19) becomes

$\displaystyle
H_{1}^{n}=\sum_{\\{k:\kappa_{\nu_{k}}=\kappa_{1}\\}}\kappa_{1}(\nu_{k}-\nu_{k-1})\sum_{j}G_{\kappa_{1}}^{ij}B_{\nu_{k}}({T^{n}_{j}})$
(21) $\displaystyle
H_{2}^{n}=\sum_{\\{k:\kappa_{\nu_{k}}=\kappa_{2}\\}}\kappa_{2}(\nu_{k}-\nu_{k-1})\sum_{j}G_{\kappa_{2}}^{ij}B_{\nu_{k}}({T^{n}_{j}})$
$\displaystyle\int_{0}^{\infty}\kappa
J^{n+1}_{\nu,i}=s_{1}\kappa_{1}S^{E}_{1}+s_{2}\kappa_{2}S^{E}_{2}+H_{1}^{n}+H_{2}^{n}$

where $s_{j}$ is the measure of the $\nu$-set where $\kappa_{\nu}=\kappa_{j}$,
$j=1,2.$

There are some limitations however: $\kappa({\bf x})$ must be a sum of
products of ${\bf x}$-functions by $\nu$\- functions:

$\kappa_{\nu}(x)=\sum_{j}\rho_{j}({\bf x})\kappa^{j}_{\nu},\hbox{ where
$\kappa^{j}_{\nu}$ takes only $N_{k}^{j}$ different values
}\\{\kappa^{j}_{k}\\}_{k}$ (22)

Let $s^{j}_{k}:=\\{\nu~{}:~{}\kappa^{j}_{\nu}=\kappa^{j}_{k}\\}.$. Then

$\displaystyle\int_{0}^{\infty}\kappa_{\nu}({\bf x}){\cal J}[S(\nu,{\bf
x})]\mathrm{d}\nu=$ (23)
$\displaystyle\frac{1}{4\pi}\sum_{i,j,l,k}\rho_{j}({\bf
x})\kappa^{j}_{k}\kappa^{i}_{l}\int_{\Omega}\left[\int_{s^{i}_{k}}S(\nu,{\bf
x}^{\prime})\mathrm{d}\nu\right]\rho_{i}(\cdot)\frac{{\rm
e}^{-\kappa^{i}_{l}\int_{[{\bf x},{\bf x}^{\prime}]}\rho_{i}(\cdot)}}{|{\bf
x}-{\bf x}^{\prime}|^{2}}\mathrm{d}x^{\prime}.$

Hence only the ${\mathcal{H}}$-matrices of kernel $\rho_{i}({\bf
x}^{\prime})\kappa^{i}_{l}\exp({-\kappa^{i}_{l}\int_{[{\bf x},{\bf
x}^{\prime}]}\rho_{i}(\cdot)})/|{\bf x}-{\bf x}^{\prime}|^{2}$ are needed.

The same decomposition holds for $S_{\nu}^{E}$.

$\displaystyle\int_{0}^{\infty}\kappa_{\nu}({\bf x})S_{\nu}^{E}({\bf
x})\mathrm{d}\nu=$ $\displaystyle\tfrac{1}{4\pi}\sum_{i,j,k,l}\rho_{i}({\bf
x})\kappa^{i}_{l}\int_{\partial\Omega}\left[\int_{s^{j}_{k}}Q_{\nu}\mathrm{d}\nu\right]{\rm
e}^{-\kappa^{j}_{k}\int_{[{\bf x}-{\bf
x}^{\prime}]}\rho(\cdot)}\left(\frac{[({\bf x}^{\prime}-{\bf x})\cdot{\bf
n}(x^{\prime})]^{+}}{|{\bf x}^{\prime}-{\bf
x}|^{2}}\right)^{2}\mathrm{d}\Gamma(x^{\prime}),$

which, after discretization, involves only the sum of products of
${\mathcal{H}}$-matrices of kernel on $\partial\Omega$: $\displaystyle{\rm
e}^{-\kappa^{j}_{k}\int_{[{\bf x}-{\bf x}^{\prime}]}\rho_{j}(\cdot)}([({\bf
x}^{\prime}-{\bf x})\cdot{\bf n}(x^{\prime})]_{-})^{2}/|{\bf x}^{\prime}-{\bf
x}|^{4}.$

###  3.2 Computations of Exponentials on a Background Grid

To compute ${\rm e}^{-\kappa^{j}_{k}\int_{[{\bf x}-{\bf
x}^{\prime}]}\rho_{j}(\cdot)}$ we use a quadrature rule for the integral:

$\int_{[{\bf x}-{\bf
x}^{\prime}]}\rho_{j}(\cdot)\approx\sum_{i=0}^{I}\rho_{j}({\bf x}_{i})|{\bf
x}_{i+1}-{\bf x}_{i}|,\quad{\bf x}_{0}={\bf x},~{}~{}{\bf x}_{I+1}={\bf
x}^{\prime}.$

To speed-up the computation of $\rho({\bf x}_{i})$, $\rho$ is interpolated on
a fine 3D Cartesian grid during the initialization phase of the computer
program.

##  4 Validation

###  4.1 Computing Resources

To assert the precision and computing time we used several meshes. The grey
case corresponds to $\kappa=0.5(1-x/2)$. The Gemini case is when
$\kappa_{\nu}=(1-x/2)\tilde{\kappa}_{\nu}$ with $\tilde{\kappa}_{\nu}$ read
from the Gemini web site. Results are on table 1 for the Chamonix valley and
on figure 7 for the flat ground case.

Table 1: Computing Times in seconds (Apple M1 8 proc) nb vertices | 4053 | 30855 | 231 796
---|---|---|---
Grey | 3 | 19 | 170
Gemini | 36 | 239 | 2027

The code has been ported on a supercomputer. It scales perfectly for the grey
case and takes 90” for 300K vertices. On the Gemini case the Newton iterations
are not yet parallelized but the rest takes 360”. The memory used is governed
by the compressed matrices. For the ${\mathcal{H}}$-matrices the compression
ratio is 29 for the volumic kernel and 9 for the surfacic kernel.

Figure 3: Compressed blocks in the ${\mathcal{H}}$-matrices, surfacic on the
left, volumic on the right. Compression is shown by colors (lighter is more
compressed, red is not compressed), numbers are the ranks of the
approximations. Neighboring blocks may not have similar ranks because they
maybe correspond to vertices far from each other in the physical space.

###  4.2 Length Scale, Values for $\kappa_{\nu}$

The Gemini measurements for $\kappa$ are used in a domain $\Omega={\bf
x}\in(0,1)\times(-0.2,3.32)\times(-3.35,0.16)$ in 10km units. The Nimbus-4
measurements [13] show that some infra-red radiations like
$\nu_{12}=\frac{3}{12}\cdot 10^{14}s^{-1}$, cross the 12km thick earth
atmosphere almost unaffected, which implies that $\kappa_{\nu_{12}}\sim 0$.

Other radiations like $\nu_{15}=0.2\cdot 10^{14}s^{-1}$ are damped by 28%
(i.e. $I_{\nu_{15}}|_{x=1}=0.28I_{\nu_{15}}|_{x=0}$, see [13]) which implies
${\rm e}^{-1.2\kappa_{\nu_{15}}}=0.28$, i.e. $\kappa_{\nu_{15}}=1.05$. This is
associated with the presence of $\texttt{CO}_{2}$ in air.

Similarly, corresponding to the presence of water in air, $I_{\nu_{20}}$ is
damped by 20%, meaning that $\kappa_{\nu_{20}}\sim 0.3$.

Therefore the Gemini data can be used without scaling. For a grey model,
$\kappa=0.5$ corresponds to the average loss from the theoretical black body
radiation at ground level to the actual measurements in space. Surely this
value is too large in the visible range, but the correlation between the
infra-red and the visible ranges is weak.

###  4.3 Temperatures on a Flat Land Exposed to Sunrays

Our purpose here is to validate the results against the stratified numerical
solutions of [8]. We begin with a grey case, $\kappa=0.5$. The domain is
$\Omega=(0,1\times(-L,L)\times(-L,L)$ above the Earth surface $S:~{}x=0$. The
radiations of the Sun cross the atmosphere unaffected; $30\%$ is reflected and
$70\%$ is absorbed and then re-emitted as a black body in all directions
(Lambertian reflection) with intensity given by(24). Since the downward travel
of the radiations is ignored, the source of light is at $x=0$, with

$Q_{\nu}({\bf x},{\bm{\omega}})=Q^{0}B_{\nu}(T_{s}){\bf
n}\cdot{\bm{\omega}},\quad Q^{0}=2\cdot 10^{-5},\quad T_{s}=1.02,~{}{\bf
n}=(1,0,0)^{T}$ (24)

The results depend on the mesh size $h$ and $L$. By increasing $L$ and the
number of vertices $n^{3}$, convergence to the stratified case is reached
(figure 5). The convergence rate is shown on figure 7 to be approximately
$h^{1.5}$. The stratified solution is computed independently by a method
sketched in Paragraph 2.0.1 and detailed in [FGOP2].

The CPU times are plotted on figure 7 for 5 meshes with N= 616, 4056, 30855,
80898, 231796 vertices respectively. It confirms the $N\log N$ growth.

A similar exercise is done with $\kappa_{\nu}=0.5(1-x/2)$. Comparison with the
stratified case is possible after a change of variable in the stratified code.
Results are shown on figure 5.

.

$0$$0.2$$0.4$$0.6$$0.8$$1$$-20$$-10$$0$$10$$20$Altitude kmTemperature CL=3,
n=5L=5, n=5L=10, n=7Stratified

Figure 4:  Temperatures $T(x,0,0)$ versus altitude $x$ when
$\kappa_{\nu}=\tfrac{1}{2}$ for some values of $L$ and $n$. Convergence to the
stratified 1D solution is visible.

.

$0$$0.2$$0.4$$0.6$$0.8$$1$$-20$$-10$$0$$10$$20$$30$Altitude kmTemperature
CL=10, n=15L=10,n=5L=10, n=10Stratified

Figure 5:  Temperatures when $\rho\tilde{\kappa}=\tfrac{1}{2}(1-0.5x)$:
comparison with the stratified solution.

.

$0$$0.5$$1$$1.5$$2$$2.5$$\cdot 10^{5}$$0$$200$$400$$600$$800$number of
verticesCPU time (seconds)CPU$y=1.5x$

Figure 6:  CPU time to compute the grey case on a flat ground, indicating an
$N\ln N$ growth where $N$ is the number of vertices in the mesh.

.

$1.6$$1.8$$2$$2.2$$2.4$$2.6$$2.8$$1.5$$2$$2.5$$3$$3.5$$\log
n$$-\log\frac{\|T-T_{strat}\|}{\|T_{strat}\|}$Relative log error$y\sim 1.5x$

Figure 7:  Log-log plot of the discretization error versus mesh size; the
“exact” solution is the stratified one computed with many points.

###  4.4 The Chamonix Valley: the Grey Case

The domain is a portion of the atmosphere
$\Omega=[h(y,z),H]\times[y_{m},y_{M}][z_{m},z_{M}]$ above the Earth surface
$S:~{}x=h(y,z)$. The radiations of the Sun cross the atmosphere unaffected;
$30\%$ is reflected and $70\%$ is absorbed and re-emitted with intensity (24),
except that ${\bf n}$ is the normal of $S$. We can account for the hour of the
day by rotating accordingly all normals.

$Q_{\nu}([h(y,z),y,z]^{T},{\bm{\omega}})=Q_{0}({\bf x})B_{\nu}(T_{s}){\bf
n}\cdot{\bm{\omega}},\quad T_{s}=1.02,$

We have simulated the case $Q_{0}=2\cdot 10^{-5}$ (no snow) and the case with
snow above 2500m:

$Q_{0}({\bf x})=0.5\left(\beta+(1-\beta){\bf
1}_{x<h_{snow}}\right),~{}~{}\beta=0.3,~{}~{}h_{snow}=0.25,~{}~{}{\bf x}\in
S.$

It means that in case of snow the light source is 30% of what it is without
snow, the rest was reflected in the visible range before re-emission. Note
that glaciers below 2500m are ignored.

As above we begin with a $\kappa_{\nu}=0.5(1-x/2)$ which takes into account
the rarefaction of air. The mesh is shown on figure 8.

Temperatures and radiations are computed in the morning and evening when the
sun is at $\pm 45^{o}$ from the vertical and inclined towards the South by
$20^{0}$. Results are shown on figure 8. On all such figures the altitude has
been multiplied by 2 to enhance the graphics. It can be seen that the
temperature is high on the slopes exposed to the Sun. It is seen also that the
temperature is much cooler where there is snow. Temperatures are clearly too
low on high mountains and perhaps too hot also on sunlit slopes. Still it is
reasonable for a grey model. The Mont-Blanc is on the top left. Recall that
points near the border of the domain receive only half the light because of
clipping, so it is unrealistically cooler near the border.

To check convergence we compared a lower solution with a higher one. To
generate a lower solution initial reduced temperature in the algorithm was set
to 0.01. To generate an upper solution, $T^{0}$ was set to 0.12. Instead of 7
iterations, 15 were needed to reduce $|T^{n+1}-T^{n}|$ below $10^{-11}$. Then,
no visible difference could be seen between the two results.

Figure 8:  All 3 computations are with $\kappa_{\nu}=0.5(1-x/2)$ and a clear
sky (no cloud). Top: Reduced temperatures at noon in a vertical cut of the
domain (the mesh has 36K vertices). Middle: Ground level temperatures in the
morning: the Sun is at $45^{o}$ from vertical in (x,y) plane and $-20^{0}$ in
the (x,z) plane. (Est is the lower side of the square, South is the left
side). Bottom :Temperatures in the evening ((sun at $-45^{o}$ from vertical in
(x,y) plane and $-20^{0}$ in the (x,z) plane; West is the upper side of the
square). Both meshes have 360K vertices.

####  4.4.1 Influence of the snow

The snow is very important for the temperature distribution at the ground
level. When there is no snow temperatures at the ground level in the Chamonix
region at noon is $2^{o}C$ hotter in the valley and very much hotter in the
mountain, above zero (figure 9).

Figure 9: Influence of the snow and the mesh: Ground temperatures at noon
without snow. It is $2^{o}C$ hotter in the valley and very much hotter in the
mountain, naturally. On top, the mesh has 36K vertices, and below, it has ten
times more vertices. the two pictures differ mostly near the border of the
domain where the computations are wrong anyway because the points receive only
one half of the sunshine other points receive (and one fourth at the corners).

####  4.4.2 Clouds

Figure 10: Ground temperature in the Chamonix region at noon without (top) and
with (bottom) the cloud. It is hotter in the cloud but colder in the valley.

Let the cloud be a cylinder centered in the middle of the domain, between
$x=2000$m and $x=8000$m. There, $\kappa_{\nu}({\bf x})$ is multiplied by
$1.5$:

$\kappa_{\nu}=0.5(1-0.5x)(1+0.5{\bf
1}[x\in(0.2,0.8):(y-1.5)^{2}+(z+1.5)^{2}<0.5]);$

On Figure 10 ground temperatures are shown with and without the cloud in a
zoom region with the finer mesh. It is hotter on the mountain inside the cloud
but colder in the valley.

Figure 13 shows the temperature functions of altitude with and without the
cloud (with snow) with variable density of air.

###  4.5 The Chamonix Valley in the Non-Grey Case

As in [8] $\nu\mapsto\kappa_{\nu}$ is read from the Gemini measurements web
site.

To be sure that $\int_{0}^{\infty}B_{\nu}(T)=\sigma T^{4}$ with a good
precision it is necessary to extend $\tilde{\kappa}_{\nu}$ for
$\nu\in(0.01,0.3$; we set it to be equal to the last available Gemini point.

Then, as shown on Figure 2, $\nu\mapsto\tilde{\kappa}_{\nu}$ is approximated
by the nearest step function that takes only 10 values:
$\texttt{round}(10\kappa)/10$.

As explained above, by this trick we need to compute only $2\times 10$
${\mathcal{H}}$-matrices, even though $\nu$-integrals are computed with 683
quadrature points.

Temperatures at noon with snow and no cloud were computed with the Gemini data
of figure 2. Results are shown on the top picture of figure 11.

Figure 11: Ground temperatures at noon with snow (no cloud) with $\kappa$
given by the Gemini data. Below, same situation, but the $\kappa_{\nu}$ from
the Gemini data is put to $1-x/2$ in the range $\nu\in(c/18,c/14)$ to simulate
the absorption by $\texttt{CO}_{2}$. Bottom: same but with
$\kappa_{\nu}=1-x/2$, $\nu\in(c/3,c/1.5)$ to simulate the absorption by
$\texttt{CH}_{4}$. The color map is the same for all 3. Bluer is colder,
browner is hotter.

Then these data were modified in the range $\nu\in(c/18,c/14)$ where
$\kappa_{\nu}$ is set to 1. It simulates grossly the effect of
$\texttt{CO}_{2}$ which renders the atmosphere opaque in these frequencies.
The results are displayed on the middle picture of figure of 11 and on figures
13 and 13. Note that the temperatures are colder.

A similar exercise is made but with $\kappa_{\nu}$ set to 1 in the range
$\nu\in(c/3,c/1.5)$ to simulate an effect in the direction of
$\texttt{CH}_{4}$. Results are in the bottom picture of figure 11 and also on
figures 13 and 13. Note that the temperatures are hotter.

Hence blocking an infrared subrange can either make the day hotter or colder,
depending on the position of the subrange, which in these numerical
experiments we have childishly called $\texttt{CO}_{2}$ and $\texttt{CH}_{4}$.
The same observation was made in [8].

On figure 13 the radiation intensities are shown versus wavelength ($c/\nu$)
above Chamonix at altitude 3000m.

Note that the computing time is roughly ten time the one with a constant
$\kappa$ because we use a Lebesgue discretization with ten levels.

$0.2$$0.4$$0.6$$0.8$$1$$-40$$-20$$0$$20$$40$Altitude kmTemperature
${}^{o}C$greywith cloud greyGeminiGemini+CO2Gemini+CH4 Figure 12:
Temperatures versus altitude in the center of the domain computed first with a
fixed $\tilde{\kappa}$ with and without the cloud. Then with a variable
$\tilde{\kappa}$ (Gemini measurements) with $\texttt{CO}_{2}$ or
$\texttt{CH}_{4}$ corrections.

.

$0$$5$$10$$15$$20$$25$$30$$0$$5$$10$$15$wave length $\mu$mlight-
intensity$I_{\nu}$ Gemini$I_{\nu}$ Gemini+CO2$I_{\nu}$ Gemini+CH4

Figure 13:  light intensity versus wave length at altitude 3000m in the center
of the domain, with Gemini-$\kappa_{\nu}$ without and with $\texttt{CO}_{2}$
or $\texttt{CH}_{4}$ corrections.

###  4.6 Effect of Thermal Diffusion and Wind Velocity

####  4.6.1 Wind velocity

Convective winds due to temperature differences and modelled by a Boussinesq
approximation seem out of simulation reach in areas of several kilometers. On
the other hand high atmosphere winds $u_{0}\sim 360km/h$ are numerically
tractable if we assume that the effective viscosity is dominated by
turbulence, leading to a Reynolds number around 500. Thus the stationary
Navier-Stokes equations with fixed varying density are solved in $\Omega$ for
the scaled velocity ${\bf u}/u_{0}$ and pressure:

${\bf u}\cdot\nabla{\bf u}+\nabla p-\nu\Delta{\bf u}=0,~{}~{}\nabla\cdot{\bf
u}=0,\quad{\bf u}(H,y,z)=[0,1,0]^{T}\quad-$ (25)

and Neumann conditions on the remaining boundaries $-\nu\partial_{\bf n}{\bf
u}+p{\bf n}=0$.

A Newton method is applied to linearize the system and the Taylor-Hood Finite
Element method is used for discretization. The linear systems are solved in
parallel with the MUMPS library. Less than a dozen iterations are sufficient
for convergence in some 600 sec on the M1 (see figure 14)

Figure 14: Solution of the Navier-Stokes eqs at Re=500. Top: ${\bf u}_{1}$,
Bottom: ${\bf u}_{2}$

####  4.6.2 Simulation of the Temperature Equation

Recall that we have neglected the variations of $\rho$ in the diffusion term:

${\bf u}\cdot\nabla T-\kappa_{T}\Delta
T=A\int_{0}^{\infty}\kappa_{\nu}(J_{\nu}-B_{\nu}(T))\mathrm{d}\nu\hbox{ in
}\Omega,\quad\frac{\partial T}{\partial n}|_{\partial\Omega}=0.$

In the computations, the unit length is L=10km. So $\kappa_{T}=2\cdot
10^{-11}$[L]2/s and $\frac{A}{\bar{\rho}}=2.70\cdot 10^{7}$[L]/s. In these
units $u_{0}=0.01$[L]/s. Let us divide by $u_{0}$:

${\bf u}\cdot\nabla T-\tilde{\kappa}_{T}\Delta
T=\tilde{A}\int_{0}^{\infty}\bar{\rho}\kappa_{\nu}(J_{\nu}-B_{\nu}(T))\mathrm{d}\nu\hbox{
in }\Omega,\quad\frac{\partial T}{\partial n}|_{\partial\Omega}=0,$

where ${\bf u}$ is the solution of (25), $\tilde{\kappa}_{T}=2\cdot 10^{-9}$,
$\tilde{A}=5.4\cdot 10^{9}$, $\bar{\rho}\kappa_{\nu}\sim 0.5$.

Evidently if $T^{*}+\delta T$ is the solution where $T^{*}$ is the solution
with ${\bf u}=0$ and $\tilde{\kappa}_{T}=0$, then $\delta T$ will be very
small. This allows for linearization:

${\bf u}\cdot\nabla\delta T-\tilde{\kappa}_{T}\Delta\delta T+\delta
T\tilde{A}\int_{0}^{\infty}\bar{\rho}\kappa_{\nu}B^{\prime}_{\nu}(T^{*})\mathrm{d}\nu=\tilde{\kappa}_{T}\Delta\delta
T^{*}\,.$

The PDE is numerically out of reach with the physical values of the parameter.
Yet to validate the concept and indicate the trend of the effect of wind and
heat diffusion, we have solved it with $\tilde{\kappa}_{T}=2\cdot 10^{-3}$,
$\bar{\rho}\kappa_{\nu}=(1-x/2)/2$ and $\tilde{A}:=A^{\prime}=5.4\cdot
10^{2}$. Notice that if $\kappa_{\nu}$ does not depend on $\nu$,
$A^{\prime}\int_{0}^{\infty}\bar{\rho}\kappa_{\nu}B^{\prime}_{\nu}(T^{*})=4\sigma{T^{*}}^{3}\bar{\rho}\kappa\sim
1.55$. The results are shown on figures 16, 16.

$0.2$$0.4$$0.6$$0.8$$1$$-20$$0$$20$Altitude kmTemperature ${}^{o}C$no windwind
Figure 15:  Temperatures versus altitude above Chamonix computed with
$\kappa=0.5(1-x/2)$, with snow and heat diffusion with and without wind.

.

Figure 16:  Temperature correction $\delta T$ on the ground due to heat
diffusion, $\kappa_{\nu}=0.5(1-0.5x)$, snow and a small wind.

##  5 Conclusion

The numerical study validates the claim that heat transfer with radiation can
be solved numerically in minutes in 3D on a laptop with O(100K) mesh nodes and
a fine resolution of the discontinuities of the absorption parameter in space
and frequencies. On a massively parallel supercomputer with O(1M) mesh nodes
it can be handled in minutes too.

The method is based on an integral formula for the radiation intensity
averaged on the unit sphere of directions. The problem is reduced to an
integral equation coupled with the temperature equation with radiation as
source. An iterative scheme can be used which is mathematically shown to be
convergent and monotone. The numerical errors can be estimated from the
difference between the lower and the upper solutions.

A drastic numerical speed-up is obtained when the convolutions in the
integrals are replaced by vector products with compressed
${\mathcal{H}}$-matrices. Furthermore a small number of matrices are needed
only when the integrals are computed as Lebesgue integrals.

The method was tested on a portion of the atmosphere around the city of
Chamonix where high mountains require a full 3D unstructured mesh. The code is
precise enough to assert the differenes due to changes in $\kappa_{\nu}$ is a
subrange of frequenceies. All results seem physically reasonable but we make
no climate claim based on the numerical results. Yet we hope to have convinced
some to try the code.

##  6 Acknowledgement

Some computations were made on the machine Joliot-Curie of the national
computing center TGCC-GENCI under allocation A0120607330.

The computer code will soon be on github at

https://github.com/FreeFem/FreeFem-
sources/tree/develop/examples/mpi/chamonix.edp

## References

  * [1] C. Bardos and O Pironneau.  Radiative Transfer for the Greenhouse Effect.  submitted to SeMA J. Springer, 2021.
  * [2] Craig F. Bohren and Eugene E. Clothiaux.  Fundamentals of Atmospheric Radiation.  Wiley-VCH, 2006.
  * [3] R. Cahalan, L. Oreopoulos, A. Marshak, K. Evans, A. Davis, and et al.  Bringing together the most advanced radiative transfer tools for cloudy atmospheres.  American Meteorological Society, 86(9):1275–1294, 2005.
  * [4] S. Chandrasekhar.  Radiative Transfer.  Clarendon Press, Oxford, 1950.
  * [5] Andrew Fowler.  Mathematical geoscience.  Berlin: Springer, 2011.
  * [6] F. Golse and O. Pironneau.  Radiative transfer in a fluid.  RACSAM, Springer, dedicated to I. Diaz(To appear), 2021.
  * [7] F. Golse and O. Pironneau.  Stratified radiative transfer for multidimensional fluids.  Compte-Rendus de l’Académie des Sciences (Mécamique), to appear 2022\.
  * [8] F. Golse and O. Pironneau.  Stratified radiative transfer in a fluid and numerical applications to earth science.  SIAM J. numerical analysis, To appear, 2022.
  * [9] Richard M. Goody and Y. L. Yung.  Atmospheric Radiation. 2nd ed. Oxford: Clarendon Press, 1996.
  * [10] F Hecht.  New developments in freefem++.  J. Numer. Math., 20:251–265, 2012.
  * [11] S.D. Lord.  Earth atmosphere transmittance measurements.  Technical report, NASA Technical Memorandum 103957,, www.gemini.edu/observing/telescopes-and-sites/sites#Transmission, 1992.
  * [12] A. Marshak and A.B. Davis, editors.  3D Radiative Transfer in Cloudy Atmospheres, volume 5117 of Physics of Earth and Space Environments.  Springer, 2005.
  * [13] F. Ortenberg.  Ozone: Space vision.  Technical Report nov, Technion, Haifa, 2002.
  * [14] O. Pironneau.  A fast and accurate numerical method for radiative transfer in the atmosphere.  Compte Rendus de l’académie des sciences (Math), December 2021.
  * [15] M. Porzio and O. Lopez-Pouso.  Application of accretive operators theory to evolutive combined conduction, convection and radiation.  Rev. Mat. Iberoamericana, 20:257–275, 2004.
  * [16] S.Boerm, L. Grasedyck, and W. Hackbusch.  Hybrid cross approximation of integral operators.  Numerische Mathematik, 10(12):221–249, 2005.
  * [17] W. Zdunkowski and T. Trautmann.  Radiation in the Atmosphere. A Course in Theoretical Meteorology.  Cambridge Univ. Press, 2007.


#  Developing Autonomous Robot-Mediated Behavior Coaching Sessions with Haru

Matouš Jelínek<EMAIL_ADDRESS>0000-0002-5303-9238 University of Southern
DenmarkSønderborgDenmark6400 ,  Eric Nichols<EMAIL_ADDRESS>0000-0003-0734-6621 Honda Research Institute Japan8-1 HonchoWakoJapan and
Randy Gomez<EMAIL_ADDRESS>0000-0002-3191-6818 Honda Research
Institute JapanWakoJapan

(2024)

###### Abstract.

This study presents an empirical investigation into the design and impact of
autonomous dialogues in human-robot interaction for behavior change coaching.
We focus on the use of Haru, a tabletop social robot, and explore the
implementation of the Tiny Habits method (fogg_tiny_2019, ) for fostering
positive behavior change. The core of our study lies in developing a fully
autonomous dialogue system that maximizes Haru’s emotional expressiveness and
unique personality. Our methodology involved iterative design and extensive
testing of the dialogue system, ensuring it effectively embodied the
principles of the Tiny Habits method while also incorporating strategies for
trust-raising and trust-dampening. The effectiveness of the final version of
the dialogue was evaluated in an experimental study with human participants
(N=12). The results indicated a significant improvement in perceptions of
Haru’s liveliness, interactivity, and neutrality. Additionally, our study
contributes to the broader understanding of dialogue design in social
robotics, offering practical insights for future developments in the field.

Autonomous Robots, Behavior Change Coaching, Human-Robot Interaction,
Emotional Expressiveness, Trust Calibration

††journalyear: 2024††copyright: rightsretained††conference: Companion of the
2024 ACM/IEEE International Conference on Human-Robot Interaction; March
11–14, 2024; Boulder, CO, USA††booktitle: Companion of the 2024 ACM/IEEE
International Conference on Human-Robot Interaction (HRI ’24 Companion), March
11–14, 2024, Boulder, CO, USA††doi: 10.1145/3610978.3640583††isbn:
979-8-4007-0323-2/24/03††ccs: Human-centered computing Empirical studies in
HCI

##  1\. Introduction

The use of social robots in behavior change coaching presents a unique
intersection of human-robot interaction (HRI), psychology, and AI. Social
robots have been previously utilized in behavior change coaching, with various
studies underscoring their potential in this field (bodala_teleoperated_2021,
; arendsen_use_2010, ; jelinek_role_2021, ). Social robots offer new ways of
interaction, learning, and behavior adaptability that can complement
traditional human-led coaching methods.

Trust emerges as a fundamental element in the realm of HRI, particularly
within the scope of behavior change coaching. The concept of trust
calibration, which involves aligning a human’s trust level with a robot’s
actual capabilities (lee_trust_2004, ; hancock_meta-analysis_2011, ), plays a
pivotal role. Proper calibration of trust is crucial for accurately reflecting
a robot’s competence and reliability, thereby preventing over-reliance or
under-utilization of robotic systems. This calibration is integral for
effective and safe interactions (de_visser_towards_2020, ).

Previous work suggests that one way of raising and maintaining the level of
trust in HRI is the robot’s emotional expressions, encompassing aspects such
as emotional routines, voice modulation, and timing. These elements
significantly influence the robot’s perceived trustworthiness and level of
engagement (breazeal_emotion_2003, ; tielman_adaptive_2014, ; nass_does_2001,
).

In behavior change coaching, where personal and emotional interactions are
paramount (fogg_tiny_2019, ), the rich expressivity of a robot can profoundly
impact the coaching effectiveness. The importance of emotional expressions in
enhancing the coaching experience was demonstrated, for example in
(jelinek_role_2021-1, ), utilizing the Tiny Habit method and the social robot
Haru.

The study revealed that emotional expressions in Haru significantly enhanced
the perceived quality of the lesson and the retention of the habit practiced.
Additionally, Haru’s presence, particularly when using emotional behaviors,
led to higher confidence in participants regarding their behavior change and
increased the likelihood of habit retention compared to control groups
presented with the same information via web content (jelinek_role_2021, ).

However, this study had various limitations: they relied on a wizard operator
and used only virtual simulation of the robot. These constraints potentially
limited the depth of interaction and the naturalness of the robot’s behavior.

Our study addresses these limitations by taking a novel approach. We developed
a new fully-autonomous dialogue for Haru, aiming to harness the full emotional
potential of the physical robot, allowing for more organic and realistic
interactions, critical in the context of behavior change coaching. By
integrating emotional routines and new voice design (Nichols_voice, ) into
Haru’s autonomous dialogue, we aim to deepen the level of engagement and trust
in HRI.

Moreover, our research extends the understanding of how emotionally
intelligent interactions influence users’ perception of the robot. In doing
so, we seek to bridge the gap between the theoretical potential of social
robots in behavior change coaching and their practical application in real-
world scenarios.

##  2\. Method

In developing the new dialogue for behavior change coaching, we embraced the
Design Thinking method (noauthor_what_nodate, ). This method is renowned for
its hands-on, user-centric approach to problem-solving.

Our primary goal was to craft a dialogue that leverages the Tiny Habits
method, explicitly tailored for Haru, a social robot. We focused on utilizing
Haru’s unique capabilities to enhance the effectiveness of the dialogue in
behavior change coaching. This section will delve into the development journey
of the dialogue and highlight the key design decisions made throughout this
process.

###  2.1. Behavior change coaching

We selected The Tiny Habit Method (fogg_tiny_2019, ) as a framework for
behavior change coaching. The Tiny Habits method is grounded in several
theoretical principles from behavior change science - e.g., (ryan2000self, ;
ajzen1991theory, ; prochaska1982transtheoretical, ). It is based on Fogg’s
Behavioral model (fogg_behavior_2009, ), which states that behavior is a
product of three elements: motivation, ability, and prompt. For behavior to
occur, a person must be i) sufficiently motivated, ii) has the ability to
perform the behavior, and iii) have the proper trigger has to be present. The
Tiny Habits method focuses on addressing all those aspects. It involves
anchoring new, small habits to existing routines, simplifying aspirations into
manageable actions, and rewarding completion with positive emotions to
integrate these habits into daily life.

###  2.2. Social Robot

We utilized Haru, an experimental tabletop social robot, as our target robot.
Haru’s design, inspired by animated characters (gomez_haru_2018, ), features
two eyes with LCD screens and moving rims equipped with addressable LED strips
alongside an LED matrix mouth. Haru boasts seven degrees of freedom, including
eye tilt, rotation, inner eye movement, base rotation, and body leaning,
allowing it to create dynamic emotional signals (gomez_haru_2018, ;
gomez_haru_emotion_2020, ). Haru’s design emphasizes the upper and middle body
for expressivity, based on studies showing these areas are crucial in
conveying emotions (gomez_haru_emotion_2020, ). Without hands, Haru’s
emotional expressions are conveyed through eye and body movements, similar to
the emotional affordances of hands and arms.

Haru’s expressivity is enhanced by its custom Text-to-Speech (TTS) voice,
featuring a range of vocal genres for emotive communication. The voice
includes seven vocal genres: default, cheeky, high-energy, question, sad,
serious, and whiny, for a broad emotional range (Nichols_voice, ). Haru’s
library of over 100 routines also supports multimodal expressions, including
body and eye motions, animations, and sounds, covering a full spectrum of
emotions (gomez_haru_emotion_2020, ).

The robot is using Google Speech-to-Text API for automatic speech recognition
(ASR), and Intent Classification and Entity Recognition Models developed by
Honda Research Institute. The dialogue structure is then loaded into behavior
trees, with all the components interacting with the robot hardware using ROS
(nichols_design_2021, ). This system, which combines off-the-shelf components
with custom robot elements, was previously tested in other applications
(nichols2022hey, ), such as small talk scenarios and hospital dialogues.

###  2.3. Implementation Process

Our development, guided by Design Thinking, involved iterative refinement of
the dialogue. Initially, we assessed the existing dialogue (jelinek_role_2021,
), identifying key issues to be changed, such as the verbosity of the
dialogue, focus only on yes/no questions, or the need for a human operator to
make all the decisions in the dialogue. Our primary goals were defined as to
develop an autonomous dialogue aligned with Haru’s personality (Nichols_voice,
), incorporating trust cues and enhancing emotional expressions, which we then
ideated and converted into a first dialogue prototype. Pilot testing,
conducted over six distinct iterations with volunteers, primarily engineers
without previous knowledge of the dialogue, was instrumental in evolving the
prototype into the final dialogue. The process focused on enhancing natural
emotional responses and trust-building strategies in the dialogue flow.

We first describe the the new dialogue developed. Subsequently, we’ll delve
into specific implemented refinements, highlighting their potential
applicability in other social robotics dialogues.

###  2.4. New Dialogue Development

The final version of the dialogue, in which Haru coaches participants to adopt
a Tiny Habit, is divided into four main sections. In the first section (S1),
Haru explains the importance of finding the right aspiration for behavior
change and then provides suggestions on how to transform this aspiration into
a small daily routine - the Tiny Habit. In the second section (S2), the
participant and the robot collaboratively identify a routine that can serve as
a reminder to conduct the new habit, referred to as the anchor moment. The
third section (S3) focuses on reinforcing the newly developed Tiny Habit by
celebrating the participant’s successful completion of the new habit. In the
last section (S4), Haru summarizes the theory and offers additional help
before concluding the session.

Figure 1. Overview of the Dialogue Structure

The image is a diagram of a three-step session for defining and implementing
Tiny Habits. The first step is ”S1 - Defining Tiny Habit,” where ’Haru’
assists users in identifying a small, achievable habit tied to a larger
aspiration and offers personalized examples. The second step is ”S2 -Defining
Anchor Moment,” where ’Haru’ aids in finding a routine to serve as a reminder
for the habit. The third step is ”S3 - Defining Celebration,” where ’Haru’
helps choose a celebration to reinforce the new habit. The diagram includes
examples for each step, such as ”I will meditate for three breaths” after the
session introduction, ”After I finish my morning coffee” as an anchor moment,
and ”I will celebrate by doing a little dance” to mark the celebration. Each
step has a color-coded bar on the side indicating ”Help available” and ”Entity
capture.” The session begins with an introduction and ends with a celebration.

In the first three sections (S1-S3), the dialogue follows a consistent
structure. In the first part (P1), which users can skip, Haru presents an
example of how he addressed the material covered in the section in his own
life. In the second part (P2), participants are asked how they want to address
the section’s topic - defining the Tiny Habit, finding the right anchor
moment, or the celebration (e.g., “What’s your aspiration you’d like to work
on?”). Based on their response, the robot selects the most appropriate answer
from one of the pre-defined categories (e.g., User reply: “I want to start
running”; Haru’s response: “Awesome! A dash of exercise can spark a fitness
firework! One way to get into shape is by scaling down to doing two pushups or
putting on your running shoes to get ready.”).

In the third part (P3), the robot confirms if the subjects have understood the
concept and offers additional coaching tips. Simultaneously, the robot begins
capturing potential user responses. The last part (P4) focuses on capturing
the user’s progress (e.g., their new Tiny Habit). Haru tries to capture the
entity multiple times, offers additional guidance and examples of good
practice, and provides a default entity for adoption if participants are
unsure (e.g., “Here’s one last idea for your possible change. You seem a tad
stressed. How about focusing on relaxation? Your new Tiny Habit could be, ‘I
will take three mindful breaths.’ Does that sound like something you might use
in the future?”). If they do not confirm, the dialogue concludes
preliminarily.

If the entity is captured during the third and fourth parts of the dialogue,
Haru confirms with the user whether it was captured correctly. If the answer
is positive, the dialogue proceeds to the next section without the possibility
of returning.

###  2.5. Dialogue Refinement

####  2.5.1. Emotional and Vocal Expressivity

We leveraged Haru’s extensive library of emotional routines and a diverse
range of voice genres to create a more engaging and human-like interaction.
Our decision to implement these routines was grounded in theoretical
principles and insights from pilot testing, which revealed a preference for
routines that included sound and were shorter than 3.5 seconds, particularly
crucial in extended dialogue segments.

Haru’s emotional routines encompass the full spectrum of Ekman’s basic
emotions (ekman1984expression, ), with multiple variants for each emotional
state. We strategically assigned these routines to specific parts of the
dialogue, ensuring they did not disrupt the flow, particularly at the end of
utterances or after questions. This approach aligns with affective grounding
theory, emphasizing shared emotional understanding as vital in human-robot
interaction (jung_affective_2017, ).

Other key enhancements in our dialogue design included integration of nuanced
emotional expressions, such as back-channels like “Oh!” and nods, to indicate
understanding and engagement. Furthermore, we tried to utilize emotional
mimicry, adapting to the emotional context of interactions, informed by
research on human emotional response patterns (kraut_social_1979, ).

These features are rooted in the concept that robots can express emotions as
signals to reveal internal states, based on basic emotion theories. Our
methodology also incorporated affective storytelling, where Haru’s dialogue
reflected appropriate emotional intonations corresponding to the narrative,
further enhancing the interaction’s intuitiveness and effectiveness
(chella_emotional_nodate, ), and trustworthiness (fischer_emotion_2019, ).

####  2.5.2. Empathy

Our study emphasized empathy in the dialogue design to foster a more
sympathetic and relatable interaction. Key to this was Haru’s ability to
respond to detected problems in the interaction with comforting statements
such as “Don’t worry, I am here to help you,” “Take your time!” and “No
problem. Let’s move forward!”. This approach is rooted in theories of
empathetic design (picard_affective_1997, ; paiva_emotion_nodate, ;
maldonado_we_2005, ), highlighting the importance of empathy in effective
interaction design. By incorporating these empathetic responses, Haru could
better connect with users, enhancing the overall trust and effectiveness of
the behavior change coaching process.

####  2.5.3. Situation Awareness

The concept of situation awareness in human-computer interaction underscores
the importance of context for creating meaningful and trustworthy interactions
(lee_trust_2004, ; abowd_towards_1999, ). In social human-robot interaction, a
robot’s situation awareness can significantly affect the robot’s persuasive
abilities (fischer_sharing_nodate, ). Recognizing this, our pilot testing
revealed the potential to enhance interaction with improved situation
awareness, thus refining the dialogue flow.

Figure 2. Illustration of the Tiny Habit Dialogue

This figure shows a person interacting with the robot Haru. There is a speech
bubble above Haru with the text, ’What’s your aspiration you’d like to work
on?’ The person responds in another bubble: ’My goal is to start exercising
every day.’ Haru replies with, ’Awesome! One example of getting into shape is
to scale down to doing two pushups or putting on your running shoes to get
ready.’

  * •

Adaptive Responses and Dialogue Branching

Haru was programmed to tailor its responses based on the user’s current
environment and expressed preferences. For instance, if the user expressed a
desire to skip a story, Haru responded with, “Let’s get right to the point
then.” Tailored suggestions were also provided in various parts of the
dialogues, depending on the user’s replies. For example, when the user was
searching for the right moment for their new routine, Haru would provide
options such as “After I pour my cereal” for mornings, “After I brush my
teeth” for evenings, and “After I log into my computer” for work settings.

  * •

Branching Based on User Input

The dialogue was scripted with branches to react more effectively to user
inputs. For example, when Haru asks, “Do you like science?” and the user
replies, “Yes, it is cool,” Haru responds enthusiastically. If the reply is “I
think it is boring,” Haru adapts with a more engaging response, “Oh, so you
think it is boring? Let me tell you something.”

  * •

Personalization and Entity Capture Override

To provide personalized coaching, Haru utilized captured entities throughout
the interaction. This included using participant’s names in dialogues (e.g.,
“Oh, {name}, what a ride we had”) and incorporating them into the context of
the dialogue (e.g., “Your anchor, along with the tiny habit, might sound
something like - ’After I finish my sandwich, {Tiny Habit of the person}’”).

Recognizing the need to capture entities while also evaluating whether a
question required a yes/no response, we implemented an override system. This
allowed Haru to capture entities effectively while still enabling users to ask
for help or clarification.

###  2.6. Error Mitigation

During pilot testing, we identified the need for a robust help system. This
led to the development of a repeat function, enabling Haru to reiterate the
last dialogue section upon request (e.g., “Can you repeat that?”).
Additionally, users were given the option to seek help (e.g., “I don’t know,
can you please help?”), utilizing the dialogue system’s intent recognition
feature. Haru also acknowledges and corrects misunderstandings or errors,
making statements like, “I’m sorry. I didn’t catch that.” or “It was tricky. I
couldn’t record your anchor moment.”, using apology as a possible trust repair
strategy (esterwood_having_2022, ; jelínek2023trust, ). Such features,
including a repeat system, enhance the user experience by providing support
and fostering a sense of understanding, which is crucial in effective human-
robot interactions.

###  2.7. Procedure

To assess the effectiveness of our interventions, we conducted an experiment
at Honda Research Institute in Japan. Participants included employees and
interns proficient in English, chosen through convenience sampling. This
approach allowed us to compare results with a prior study (jelinek_role_2021,
). Participants were first provided with a written consent form, ensuring
voluntary participation and agreement with data handling. They were then
familiarized with the social robot and the behavior change session, with an
opportunity to ask questions before starting.

The sessions with a physically embodied Haru, lasting approximately 12
minutes, were recorded using dual cameras: one capturing participant reactions
and the other focusing on the robot. Post-session, participants completed a
questionnaire to evaluate the session quality. This included selected items
from the Godspeed questionnaire (bartneck_measurement_2008, ) focusing on
anthropomorphism, animacy, and likeability, alongside the MDMT trust scale
(ullman_mdmt_2020, ), and a set of open questions. We also gathered basic
demographic data and inquired about prior experiences with robots.

##  3\. Results

A group of 12 individuals with an average age of 29 years (Stdev=5.72 years)
participated in the experiment. Most participants were male, constituting 83%
of the sample, with 17% identifying as female. Regarding previous robot
experience, 25% of participants have only seen robots a few times in reality,
while the remaining 75% have had a more direct engagement, either playing or
regularly working with robots. The participants’ origins were predominantly
Asian (58%), followed by North Americans (25%) and Europeans (17%).

In data analysis, conducted using SPSS, we compared our results with the data
collected in Condition 1 of the previous experiment (jelinek_role_2021, ).
Utilizing an Independent Samples T-Test, our findings supported the hypothesis
that the changes made in our experiment had a significant effect on the
session outcomes.

In the newly developed dialogue, Haru was perceived more positively compared
to earlier versions. The statistical analysis revealed that Haru was perceived
as significantly less fake (p = 0.042), significantly more lively (p = 0.002),
and interactive (p = 0.005). These results indicate a notable shift in
participants’ perceptions, aligning with the objectives of our dialogue
enhancements.

The qualitative data from the experiment with Haru revealed a mix of positive
and negative reactions. Participants generally appreciated Haru’s expressive
characteristics, such as its positive demeanor, movements, and specific
routines like chuckling, which added to its human-like appeal. The flow of
conversation with Haru was perceived as natural, with an ability to maintain
context, and its responsiveness and timing were praised. Visual elements like
eye animations and movements also enhanced the experience. Overall, the
interaction was described as enjoyable. However, issues with the ASR system
were noted multiple times (including cases where the robot ’misunderstood’
users and followed incorrect dialogue branches), which likely impacted the
overall session evaluation.

See the experiment video for implementation examples: link.

Item | MD | p | SD
---|---|---|---
Anthropomorphism |  |  |
Fake/Neutral | -0.612 | 0.042 | 0.237
Machinelike/Humanlike | -0.100 | 0.282 | 0.288
Unconscious/Conscious | 0.067 | 0.307 | 0.315
Artificial/Lifelike | -0.017 | 0.337 | 0.310
Moving rigidly/elegantly | 0.092 | 0.267 | 0.309
Animacy |  |  |
Dead/Alive | -0.083 | 0.293 | 0.292
Stagnant/Lively | -1.117 | 0.002 | 0.244
Mechanical/Organic | -0.117 | 0.265 | 0.267
Artificial/Lifelike | -0.083 | 0.298 | 0.321
Inert/Interactive | -0.850 | 0.005 | 0.214
Apathetic/Responsive | -0.233 | 0.192 | 0.270
Likability |  |  |
Dislike/Like | -0.450 | 0.088 | 0.267
Unfriendly/Friendly | -0.250 | 0.159 | 0.230
Unkind/Kind | -0.250 | 0.139 | 0.203
Unpleasant/Pleasant | -0.267 | 0.146 | 0.226
Awful/Nice | -0.133 | 0.230 | 0.209
Incompetent/Competent | 0.405 | 0.031 | 0.208
Table 1. Summary of mean differences and statistical significance for the
measured items.

##  4\. Conclusion

Our research presents an advancement in the field of human-robot dialogue
design, especially in the context of autonomous behavior change coaching with
Haru. Through rigorous pilot testing and evaluation, we developed a refined
approach to dialogue design, enhancing situation awareness and trust-building
strategies. The implementation of a fully autonomous dialogue system, enriched
with nuanced situation awareness, robot emotiveness, and a variety of trust-
building strategies, marked a novel approach in this domain. The effectiveness
of these interventions was reflected in our statistical analysis, showcasing
the positive impacts of these innovations. However, the study also revealed
potential flaws, such as issues with the Automatic Speech Recognition (ASR)
system, which provide valuable insights for future improvements.

Our work not only contributes to the general knowledge of Haru dialogue design
but sets a precedent for future research and applications in robot-mediated
behavior coaching, highlighting the potential of empathetic, context-aware,
and trustworthy interactions in this evolving field. Furthermore, the
knowledge and insights gained from our Design Thinking (DT) approach provide a
robust foundation for designing future conversations not only for Haru,
showcasing how nuanced and effective dialogue design can significantly enhance
user experience in diverse contexts.

##  5\. Future work

For future work, we aim to utilize the developed dialogue as a dynamic
platform for continued research, investigating the specific effects of the
isolated dialogue improvements on coaching efficacy and trust. Additionally,
we might focus on long-term evaluations to understand Haru’s sustained impact
in behavior change coaching. Addressing technical challenges and further
development of the dialogue system refinement might also be future research.

##  6\. ACKNOWLEDGMENT

This project was funded by Honda Research Institute Japan and the Independent
Research Fund Denmark, grant number 1032-00311B.

## References

  * [1] Abowd, Gregory Abowd, Anind Dey, Anind K, Brown, Peter J, Davies, Nigel, Smith, Mark Sandler, Steggles, and Pete.  Towards a Better Understanding of Context and Context-Awareness.  1999\.
  * [2] Icek Ajzen.  The theory of planned behavior.  Organizational behavior and human decision processes, 50(2):179–211, 1991.
  * [3] Jeroen Arendsen, Joris Janssen, Sander Begeer, and Fred Stekelenburg.  The use of robots in social behavior tutoring for children with ASD.  2010\.  Pages: 372.
  * [4] Christoph Bartneck, Dana Kulic, Elizabeth Croft, and Susana Zoghbi.  Measurement instruments for the anthropomorphism, animacy, likeability, perceived intelligence, and perceived safety of robots.  1:71–81, 2008.
  * [5] Indu Bodala, Nikhil Churamani, and Hatice Gunes.  Teleoperated Robot Coaching for Mindfulness Training: A Longitudinal Study.  2021\.
  * [6] Cynthia Breazeal.  Emotion and sociable humanoid robots.  59(1):119–155, 2003.
  * [7] A Chella.  An emotional storyteller robot.  2007\.
  * [8] Ewart J. de Visser, Marieke M. M. Peeters, Malte F. Jung, Spencer Kohn, Tyler H. Shaw, Richard Pak, and Mark A. Neerincx.  Towards a theory of longitudinal trust calibration in human–robot teams.  12(2):459–478, 2020.
  * [9] Paul Ekman.  Expression and the nature of emotion.  Approaches to emotion, 3(19):344, 1984.
  * [10] Connor Esterwood and Lionel P Robert.  Having the right attitude: How attitude impacts trust repair in human-robot interaction.  page 10.
  * [11] Kerstin Fischer, Lars Christian Jensen, and Nadine Petersen.  Sharing the situational context influences a speaker’s (a robot’s) persuasiveness.  page 19, 2022.
  * [12] Kerstin Fischer, Malte Jung, Lars Jensen, and Maria aus der Wieschen.  Emotion Expression in HRI - When and Why. 10.1109/HRI.2019.8673078.
  * [13] B. J. Fogg.  Tiny Habits: The Small Changes That Change Everything.  Houghton Mifflin Harcourt, 2019.  Google-Books-ID: 5TmbDwAAQBAJ.
  * [14] Bj Fogg.  A behavior model for persuasive design.  In Proceedings of the 4th International Conference on Persuasive Technology, pages 1–7. ACM, 2009.
  * [15] Interaction Design Foundation.  What is design thinking?
  * [16] Randy Gomez, Deborah Szapiro, Kerl Galindo, and Keisuke Nakamura.  Haru: Hardware design of an experimental tabletop robot assistant.  In Proceedings of the 2018 ACM/IEEE International Conference on Human-Robot Interaction, pages 233–240. ACM, 2018.
  * [17] Randy Gomez, Deborah Szapiro, Luis Merino, and Keisuke Nakamura.  A holistic approach in designing tabletop robot’s expressivity.  In 2020 IEEE International Conference on Robotics and Automation (ICRA), pages 1970–1976, 2020.
  * [18] Peter A. Hancock, Deborah R. Billings, Kristin E. Schaefer, Jessie Y. C. Chen, Ewart J. de Visser, and Raja Parasuraman.  A meta-analysis of factors affecting trust in human-robot interaction.  53(5):517–527, 2011.
  * [19] Matouš Jelínek and Kerstin Fischer.  The role of a social robot in behavior change coaching.  In Companion of the 2021 ACM/IEEE International Conference on Human-Robot Interaction, pages 434–438. ACM, 2021.
  * [20] Matouš Jelínek and Kerstin Fischer.  The role of emotional expression in behavior change coaching by a social robot.  In Raian Ali, Birgit Lugrin, and Fred Charles, editors, Persuasive Technology, Lecture Notes in Computer Science, pages 193–199. Springer International Publishing, 2021.
  * [21] Matouš Jelínek and Kerstin Fischer.  Trust regulation in social robotics: From violation to repair.  CHI Conference on Human Factors in Computing Systems, SARs: TMI Workshop, 2023.
  * [22] Malte Jung.  Affective Grounding in Human-Robot Interaction.  2017\.
  * [23] Robert E. Kraut and Robert E. Johnston.  Social and emotional messages of smiling: An ethological approach.  37(9):1539–1553, 1979.  Place: US Publisher: American Psychological Association.
  * [24] John D. Lee and Katrina A. See.  Trust in automation: Designing for appropriate reliance.  46(1):50–80, 2004.  Publisher: SAGE Publications Inc.
  * [25] Heidy Maldonado, Jong-Eun Lee, Scott Brave, Clifford Nass, Hiroshi Nakajima, Ryota Yamada, Kimihiko Iwamura, and Yasunori Morishima.  We learn better together: enhancing eLearning with emotional characters. 2005\.  Pages: 417.
  * [26] Clifford Nass and Kwan Min Lee.  Does computer-synthesized speech manifest personality? experimental tests of recognition, similarity-attraction, and consistency-attraction.  7(3):171–181, 2001.  Place: US Publisher: American Psychological Association.
  * [27] Eric Nichols, Leo Gao, Yurii Vasylkiv, and Randy Gomez.  Design and analysis of a collaborative story generation game for social robots.  3:674333, 2021.
  * [28] Eric Nichols, Sarah Rose Siskind, Levko Ivanchuk, Guillermo Pérez, Waki Kamino, Selma Šabanović, and Randy Gomez.  Hey haru, let’s be friends! using the tiers of friendship to build rapport through small talk with the tabletop robot haru.  In 2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pages 6101–6108. IEEE, 2022.
  * [29] Eric Nichols, Sarah Rose Siskind, Waki Kamino, Selma Šabanović, and Randy Gomez.  Iterative design of an emotive voice for the tabletop robot haru.  In Haizhou Li, Shuzhi Sam Ge, Yan Wu, Agnieszka Wykowska, Hongsheng He, Xiaorui Liu, Dongyu Li, and Jairo Perez-Osorio, editors, Social Robotics, pages 362–374, Cham, 2021. Springer International Publishing.
  * [30] Ana Paiva, Iolanda Leite, and Tiago Ribeiro.  Emotion modelling for social robots.  2014\.
  * [31] Rosalind W. Picard.  Affective computing.  Affective computing. The MIT Press, 1997.  Pages: xii, 292.
  * [32] James O Prochaska and Carlo C DiClemente.  Transtheoretical therapy: Toward a more integrative model of change.  Psychotherapy: theory, research & practice, 19(3):276, 1982.
  * [33] Richard M Ryan and Edward L Deci.  Self-determination theory and the facilitation of intrinsic motivation, social development, and well-being.  American psychologist, 55(1):68, 2000.
  * [34] Myrthe Tielman, Mark Neerincx, John-Jules Meyer, and Rosemarijn Looije.  Adaptive emotional expression in robot-child interaction.  In Proceedings of the 2014 ACM/IEEE international conference on Human-robot interaction, HRI ’14, pages 407–414. Association for Computing Machinery, 2014.
  * [35] Daniel Ullman and Bertram F Malle.  MDMT: Multi-dimensional measure of trust.  page 3, 2020.