# When the ends don’t justify the means: Learning a treatment strategy to
prevent harmful indirect effects

Kara E. Rudolph
Department of Epidemiology  The authors gratefully acknowledge R00DA042127
from the National Institute on Drug Abuse Columbia University
and
Iván Díaz
Department of Population Health Sciences  Weill Cornell Medicine

###### Abstract

There is a growing literature on finding so-called optimal treatment rules,
which are rules by which to assign treatment to individuals based on an
individual’s characteristics, such that a desired outcome is maximized. A
related goal entails identifying individuals who are predicted to have a
harmful indirect effect (the effect of treatment on an outcome through
mediators) even in the presence of an overall beneficial effect of the
treatment on the outcome. In some cases, the likelihood of a harmful indirect
effect may outweigh a likely beneficial overall effect, and would be reason to
caution against treatment for indicated individuals. We build on both the
current mediation and optimal treatment rule literature to propose a method of
identifying a subgroup for which the treatment effect through the mediator is
harmful. Our approach is nonparametric, incorporates post-treatment variables
that may confound the mediator-outcome relationship, and does not make
restrictions on the distribution of baseline covariates, mediating variables
(considered jointly), or outcomes. We apply the proposed approach to identify
a subgroup of boys in the Moving to Opportunity housing voucher experiment who
are predicted to have harmful indirect effects, though the average total
effect is beneficial.

Keywords: dynamic treatment regime, optimal individualized treatment regime,
cross-validation, causal inference, optimal rule, mediation, interventional
indirect effect

##  1 Introduction

When individuals differ in their responses to a treatment or intervention, it
may be desirable to use individual-level information to predict for whom
treatments/interventions may work well versus not. This goal has fueled the
growing literature on finding so-called optimal individualized treatment
rules/regimes, which assign treatment based on an individual’s characteristics
such that a desired outcome is maximized (Murphy, 2003). For a single point in
time and a binary treatment, the optimal treatment rule is often defined as
the maximizer over all possible treatment rules of the counterfactual mean
outcome observed under a hypothetical implementation of the rule (Murphy,
2003; Robins, 2003).

Nabi et al. (2018) and Shpitser and Sherman (2018) recently extended these
methods to finding a treatment rule that optimizes a path-specific effect.
Nabi et al. (2018) focus, in particular, on optimizing a direct effect (an
effect of treatment on an outcome, not through a mediator); in their
application, the direct effect of interest is the effect of a treatment not
operating through adherence. A related goal entails identifying individuals
who are predicted to have a harmful indirect effect (the effect of treatment
on an outcome through mediators). Such a goal would be of interest if the
process by which an effect of a treatment or intervention might occur may be
traumatic or otherwise harmful, or if the harmful indirect effect is
indicative of more general harmful social processes, and consequently
desirable to avoid. It is possible for harmful indirect effects to be present
even if the overall effect of treatment on the outcome is beneficial,
sometimes referred to as inconsistent mediation (MacKinnon et al., 2000). In
some cases, the likelihood of a harmful indirect effect may outweigh a likely
beneficial overall effect and would be reason to caution against treatment for
indicated individuals.

The Moving to Opportunity study (MTO) presents an example where it may be of
interest to identify those with potentially harmful indirect effects. In MTO,
families were randomized to receive a Section 8 housing voucher that would
allow them to move out of public housing and into a rental on the private
market (Kling et al., 2007). Participants were followed for 10-15 years and a
broad array of outcomes related to health, education and income were assessed.
Among the children, boys in the voucher group were generally found to have
negative effects on mental health and risk behaviors (though this was not true
for all such outcomes) whereas girls were generally found to have positive
effects (Clampet-Lundquist et al., 2011; Sanbonmatsu et al., 2011; Schmidt et
al., 2017; Kessler et al., 2014; Rudolph et al., 2018; Kling et al., 2007;
Rudolph et al., 2020). Much of the negative effects were explained by
mediators related to the neighborhood and school environments and instability
of the social environment (Rudolph et al., 2020). However, measures of the
school and neighborhood environments appeared to generally improve among those
receiving the intervention, indicating that the “objective” mediator measures
of a school’s academic performance, poverty level, etc., may in fact be
surrogates for harmful processes like discrimination and increased risk of
suspensions/expulsions for some (Rudolph et al., 2020). In the context of MTO,
one could argue that a slightly beneficial long-term overall effect of housing
voucher receipt on reduced risk of a particular outcome, e.g., adolescent
alcohol use, may be outweighed by harmful indirect effects, as these may
indicate the occurrence of harmful social processes more generally.

Such scenarios motivate identifying those predicted to have harmful indirect
effects, and for that subgroup, treatment may be reconsidered. In this sense,
such subgroup identification is descriptive rather than prescriptive. Although
the approach we propose is mathematically grounded in the optimal
individualized treatment rule literature, we avoid using the phrase “treatment
rule” in describing it, because our subgroup identification is not a treatment
recommendation (in part because the overall treatment effect may still be
positive). A decision to treat (or not) individuals in this subset would
likely be based on a comprehensive assessment involving aspects like the
predicted direct effect as well as ethics, costs, and qualitative measures.

In this paper, we build on both the current mediation and optimal
individualized treatment rule literature to propose a method of identifying a
subgroup who may have an unintended harmful indirect effect, and then estimate
how the population path-specific effects would hypothetically change if the
identified subgroup were not treated. Our approach is nonparametric,
incorporates an ensemble of machine learning algorithms in model fitting,
incorporates post-treatment variables that may confound the mediator-outcome
relationship, and does not make restrictions on the distribution of baseline
covariates, mediating variables (considered jointly), or outcomes. We apply
the proposed approach to identify a subgroup of boys in the MTO housing
voucher experiment who are predicted to have harmful indirect effects, though
the average total effect on adolescent alcohol use is beneficial.

##  2 Notation

Let $O=(W,A,Z,M,Y)$ represent observed data, where $W$ is a vector of baseline
variables; $A$ is a binary treatment representing housing voucher receipt; $M$
represents mediating variables related to the school, neighborhood, and social
environments; $Y$ represents the outcome of alcohol use in adolescence, and
$Z$ represents post-treatment variables that may confound the $M\rightarrow Y$
relationship (in the case of MTO, moving with the housing voucher). Let
$O_{1},...,O_{n}$ represent the sample of $n$ i.i.d. observations of $O$. We
use capital letters to denote random variables and lowercase letters to denote
realizations of those variables. Let $\mathsf{P}$ represent the distribution
of $O$. We assume the data are generated according to the following
nonparametric structural equation model (Pearl, 2000):

$W=f_{W}(U_{W});A=f_{A}(W,U_{A});Z=f_{Z}(A,W,U_{Z});\\\
M=f_{M}(Z,A,W,U_{M});Y=f_{Y}(M,Z,A,W,U_{Y}),$ (1)

where each $U$ represents unmeasured, exogenous factors and each $f$
represents an unknown, deterministic function.

Our method to identify subsets of individuals at risk of a harmful indirect
effect are defined in terms of counterfactual outcomes under interventions on
the treatment and mediator. Specifically, we define a counterfactual outcome
under interventions that set the treatment and mediator to $(A,M)=(a,m)$ as
$Y_{a,m}=f_{Y}(m,Z_{a},a,W,U_{Y})$, where $Z_{a}=f_{Z}(a,W,U_{Z})$ is the
counterfactual variable $Z$ observed under an intervention setting $A=a$.
Direct and indirect effects are not generally point-identified in the presence
of a variable $Z$ that confounds the mediator-outcome relation and that is
affected by treatment (Avin et al., 2005). We overcome this problem using so-
called interventional effects (VanderWeele et al., 2014; Vansteelandt and
Daniel, 2017; Zheng and van der Laan, 2017; Lok, 2016, 2019; Rudolph et al.,
2017; Didelez et al., 2006), which rely on interventions on that set the
mediator to a random draw from its counterfactual distribution conditional on
baseline variables. Specifically, we let $M_{a}=f_{M}(Z_{a},a,W,U_{M})$, and
use $G_{a}$ to denote a random draw from the distribution of $M_{a}$
conditional on $W$.

We define $\mathsf{P}f=\int f(o)\mathrm{d}\mathsf{P}(o)$ for a given function
$f(o)$. We assume $\mathsf{P}$ is continuous with respect to some dominating
measure $\nu$ and let $\mathsf{p}$ denote the corresponding probability
density function. We will also use the following notation:
$\mathsf{b}(a,z,m,w)$ denotes $\mathsf{E}(Y\mid A=a,Z=z,M=m,W=w)$;
$\mathsf{g}(a\mid w)$ denotes $\mathsf{P}(A=a\mid W=w)$; $\mathsf{e}(a\mid
m,w)$ denotes $\mathsf{P}(A=a\mid M=m,W=w)$; $\mathsf{q}(z\mid a,w)$ denotes
$\mathsf{P}(Z=z\mid A=a,W=w)$; and $\mathsf{r}(z\mid a,m,w)$ denotes
$\mathsf{P}(Z=z\mid A=a,M=m,W=w)$.

##  3 Identification of the subgroup with a predicted harmful indirect effect

###  3.1 Estimand and identification

Throughout, we use the example from MTO in which a reduction in risk of
subsequent adolescent alcohol use is the outcome of interest; thus, negative
risk differences are considered “beneficial” and positive risk differences are
considered “harmful”. The definition of “harmful” for what follows should be
changed to reflect the particular research question.

As noted above, we use population interventional direct and indirect effects.
We only provide their definition here, as they have have been described
extensively elsewhere (Petersen et al., 2006; van der Laan and Petersen, 2008;
Zheng and van der Laan, 2012; VanderWeele et al., 2014; Rudolph et al., 2017;
VanderWeele and Tchetgen Tchetgen, 2017; Díaz et al., 2019). The population
interventional indirect effect,
$\mathsf{E}(Y_{1,G_{1}})$-$\mathsf{E}(Y_{1,G_{0}})$, compares the expected
average counterfactual outcomes under hypothetical interventions in which the
treatment is fixed but the mediator is changed from $G_{1}$ to $G_{0}$. This
corresponds to the effect of treatment on the outcome operating through $M$,
and is discussed further in VanderWeele and Tchetgen Tchetgen (2017). In terms
of our MTO research question, the population interventional indirect effect is
interpreted as the average difference in predicted risk of adolescent alcohol
use setting voucher to be received and drawing mediator values from their
counterfactual joint distribution if the voucher had been received versus had
it not been received; or more simply, the effect of voucher receipt on the
long-term risk of adolescent alcohol use that operates through aspects of the
school, neighborhood, and social environments.

The subgroup predicted to not have a harmful indirect effect (where “harmful”
is indicated by a positive effect, as described above) can be identified by
the target parameter

$\mathsf{d}(v)=\mathds{1}(\mathsf{B}(v)\leq 0),$

where $\mathsf{B}(v)$ is the average indirect effect conditional on a subset
of baseline variables $V\subseteq W$, and is defined as

$\mathsf{B}(v)=\mathsf{E}(Y_{1,G_{1}}-Y_{1,G_{0}}\mid V=v).$

In what follows we let $W^{c}=W\backslash V$. The function $\mathsf{B}$ is
often referred to as a “blip” function (Robins, 2004), and is a predictive
function that takes covariates $V$ as input and outputs the conditional
indirect effect of treatment on an additive scale. The function
$\mathsf{d}(v)$ is then an indicator of a predicted nonharmful indirect effect
for covariate profile $V=v$. In this case, if $\mathsf{B}(v)\leq 0,$ then
individuals in strata $V=v$ are indicated as not having a predicted harmful
effect; if $\mathsf{B}(v)>0,$ then individuals are indicated as having a
predicted harmful effect.

We introduce the following assumptions, which will allow us to identify the
causal parameters $\mathsf{B}(v)$ and $\mathsf{d}(v)$:

######  Assumption 1 (No unmeasured confounders of the $A\rightarrow Y$
relation).

$Y_{1,m}\mbox{$\perp\\!\\!\\!\perp$}A\mid W$

######  Assumption 2 (No unmeasured confounders of the $A\rightarrow M$
relation).

$M_{a}\mbox{$\perp\\!\\!\\!\perp$}A\mid W$, for $a\in\\{0,1\\}$.

######  Assumption 3 (No unmeasured confounders of the $M\rightarrow Y$
relation).

$Y_{1,m}\mbox{$\perp\\!\\!\\!\perp$}M\mid(W,A,Z)$.

######  Assumption 4 (Positivity).

$\mathsf{p}(W)>0$ implies $\mathsf{p}(a\mid W)>0$, $\mathsf{p}(M\mid
a^{\star},W)>0$ and $\mathsf{p}(Z\mid a^{\prime},W)>0$ imply $\mathsf{p}(M\mid
a^{\prime},Z,W)>0$ with probability one for $(a^{\prime},a^{\star})=(1,1)$ and
$(a^{\prime},a^{\star})=(1,0)$.

In the context of the conditional indirect effect of MTO housing voucher
receipt on subsequent risk of adolescent alcohol use through the school,
neighborhood and social environments, Assumptions 1 and 2 are expected to
hold, as the intervention is randomized. Assumption 3 may not hold, most
likely due to unmeasured post-randomization confounding variables, such as
those related to changes in parental employment, income, parent-child
dynamics, etc. However, we do include a large number of baseline covariates
related parental socioeconomic status, motivations for enrolling in MTO, and
relationships; child characteristics; and neighborhood characteristics. We
also include an indicator representing moving with the housing voucher, which
could be one such post-treatment confounder. Lastly, we include numerous
mediating variables capturing aspects of the school, neighborhood and social
environments. Any unmeasured variables would need to contribute to mediator-
outcome confounding independently of the aforementioned variables to violate
Assumption 3.

Under Assumptions 1-4, the conditional indirect effect is identified from the
observed data as follows. First, VanderWeele and Tchetgen Tchetgen (2017) show
that for any values $(a^{\prime},a^{*})\in\\{0,1\\}^{2}$,
$\mathsf{E}(Y_{a^{\prime},G_{a^{*}}}\mid V=v)$ is identified as

$\mathsf{E}(Y_{a^{\prime},G_{a^{*}}}\mid
V=v)=\int\mathsf{b}(a^{\prime},z,m,w)\mathsf{q}(z\mid
a^{\prime},w)\mathsf{p}(m\mid a^{*},w)\mathsf{p}(w^{c}\mid
v)\mathrm{d}\nu(w^{c},z,m).$

And thus we have

$\mathsf{B}(v)=\int\mathsf{b}(1,z,m,w)\mathsf{q}(z\mid 1,w)\\{\mathsf{p}(m\mid
1,w)-\mathsf{p}(m\mid 0,w)\\}\mathsf{p}(w^{c}\mid v)\mathrm{d}\nu(w^{c},z,m).$
(2)

#### Alternative interpretation in terms of the optimal dynamic treatment
rule.

Although in the context of identifying boys in MTO who are predicted to have
harmful indirect effects, we believe it would be inappropriate to interpret
$\mathsf{d}(v)$ as a treatment rule, we recognize that other contexts exist in
which such an interpretation may be justified. In these cases, $\mathsf{d}(v)$
may be alternatively defined in terms of an optimal individualized dynamic
treatment rule:

$\mathsf{d}(v)\in\operatorname*{arg\,min}_{\mathsf{d}^{\prime}\in\mathcal{D}}\mathsf{E}(Y_{\mathsf{d}^{\prime},G_{\mathsf{d}^{\prime}}})-\mathsf{E}(Y_{\mathsf{d}^{\prime},G_{0}}),$

where, for any $\mathsf{d}^{\prime}$ and $\mathsf{d}^{\star}$,
$Y_{\mathsf{d}^{\prime},G_{\mathsf{d}^{\star}}}$ is the counterfactual outcome
in a hypothetical world where treatment is set to the rule
$\mathsf{d}^{\prime}(V)\in\\{0,1\\}$ and the mediator is set to a random draw
$G_{\mathsf{d}^{\star}}$ from the distribution of $M$ conditional on
$A=\mathsf{d}^{\star}(V)$. Here $\mathcal{D}$ is the space of all functions
that map the covariates $V$ to a treatment decision rule in $\\{0,1\\}$.

###  3.2 Estimation

As pointed out in the Introduction, estimation of $\mathsf{B}(v)$ and
$\mathsf{d}(v)$ are equivalent to estimation of conditional average effects
(here, the conditional interventional indirect effect) and estimation of
optimal treatment rules. Consequently, estimation techniques such as
Q-learning (Murphy, 2003; Qian and Murphy, 2011; Moodie et al., 2012; Laber et
al., 2014), outcome-weighted learning (OWL) (Zhang et al., 2012; Zhao et al.,
2012, 2015), and doubly robust techniques (Luedtke and van der Laan, 2016;
Díaz et al., 2018; Kennedy, 2020) may be used. Q-learning is a regression-
based approach related to g-computation that relies on sequential regression
formulas. Q-learning is not directly applicable to our problem because
Equation 2 does not readily yield a sequential regression representation. OWL
uses inverse probability weights to recast the problem of estimating
$\mathsf{d}(v)$ as a weighted classification problem. Although we could use
OWL to estimate $\mathsf{B}(v)$ and $\mathsf{d}(v)$ here, we instead choose to
use a doubly robust approach that combines regressions with inverse
probability weights to obtain an estimator that remains consistent under
certain configurations of inconsistent estimation of nuisance parameters. This
approach relies on so-called unbiased transformations, which we define below.

######  Definition 1 (Unbiased transformation).

The function $D(o)$ is an unbiased transformation for $\mathsf{B}(v)$ if it
satisfies $\mathsf{E}[D(O)\mid V=v]=\mathsf{B}(v)$.

In particular, we use a multiply robust unbiased transformation (Rubin and van
der Laan, 2007). For the case of optimal treatment rules, multiply robust
unbiased transformations are constructed using the efficient influence
function (EIF) of the average treatment effect (Luedtke and van der Laan,
2016; Díaz et al., 2018; Kennedy, 2020). We take a similar approach here.
Specifically, the EIF for the counterfactual mean
$\mathsf{E}(Y_{a^{\prime},G_{a}^{*}})$ under Assumptions 1-4 was derived by
Díaz et al. (2019) as follows. Define

$\displaystyle\mathsf{h}(z,m,w)$
$\displaystyle=\frac{\mathsf{g}(a^{\prime}\mid w)}{\mathsf{g}(a^{\star}\mid
w)}\frac{\mathsf{q}(z\mid a^{\prime},w)}{\mathsf{r}(z\mid
a^{\prime},m,w)}\frac{\mathsf{e}(a^{\star}\mid m,w)}{\mathsf{e}(a^{\prime}\mid
m,w)}$ $\displaystyle\ \mathsf{u}(z,w)$
$\displaystyle=\mathsf{E}\left\\{\mathsf{b}(a^{\prime},Z,M,W)\mathsf{h}(Z,M,W),\mid\,Z=z,A=a^{\prime},W=w\right\\},$
$\displaystyle\mathsf{v}(w)$
$\displaystyle=\mathsf{E}\left\\{\int_{\mathcal{Z}}\mathsf{b}(a^{\prime},z,M,W)\mathsf{q}(z\mid
a^{\prime},W)\mathrm{d}\nu(z)\,\mid\,A=a^{\star},W=w\right\\}.$

Let
$\eta=(\mathsf{g},\mathsf{e},\mathsf{q},\mathsf{r},\mathsf{b},\mathsf{u},\mathsf{v})$
denote a vector of nuisance parameters. Then the EIF for
$\mathsf{E}(Y_{a^{\prime},G_{a}^{*}})$ under Assumptions 1-4 is equal to

$\displaystyle D_{\eta}^{(a^{\prime},a^{\star})}(o)$
$\displaystyle=\frac{\mathds{1}\\{a=a^{\prime}\\}}{\mathsf{g}(a^{\prime}\mid
w)}\mathsf{h}(z,m,w)\\{y-\mathsf{b}(a^{\prime},z,m,w)\\}$
$\displaystyle+\frac{\mathds{1}\\{a=a^{\prime}\\}}{\mathsf{g}(a^{\prime}\mid
w)}\\{\mathsf{u}(1,w)-\mathsf{u}(0,w)\\}\left\\{z-\mathsf{q}(1\mid
a^{\prime},w)\right\\}$
$\displaystyle+\frac{\mathds{1}\\{a=a^{\star}\\}}{\mathsf{g}(a^{\star}\mid
w)}\left\\{\int\mathsf{b}(a^{\prime},z,m,w)\mathsf{q}(z\mid
a^{\prime},w)\mathrm{d}\nu(z)-\mathsf{v}(w)\right\\},$
$\displaystyle+\mathsf{v}(w).$

Thus, we use $D_{\eta}(o)=D_{\eta}^{(1,1)}(o)-D_{\eta}^{(1,0)}(o)$ as an
unbiased transformation for $\mathsf{B}(v)$. For a fixed
$\eta_{1}=(\mathsf{g}_{1},\mathsf{e}_{1},\mathsf{q}_{1},\mathsf{r}_{1},\mathsf{b}_{1},\mathsf{u}_{1},\mathsf{v}_{1})$,
Díaz et al. (2019) show that $D_{\eta_{1}}(o)$ satisfies Definition 1 as long
as

  1. 1.

$\mathsf{v}_{1}=\mathsf{v}$ and either
$(\mathsf{q}_{1},\mathsf{e}_{1},\mathsf{r}_{1})=(\mathsf{q},\mathsf{e},\mathsf{r})$
or $(\mathsf{b}_{1},\mathsf{q}_{1})=(\mathsf{b},\mathsf{q})$, or

  2. 2.

$\mathsf{g}_{1}=\mathsf{g}$ and either
$(\mathsf{q}_{1},\mathsf{e}_{1},\mathsf{r}_{1})=(\mathsf{q},\mathsf{e},\mathsf{r})$
or $(\mathsf{b}_{1},\mathsf{q}_{1})=(\mathsf{b},\mathsf{q}).$

In this sense, using a multiply robust unbiased transformation provides
robustness to misspecification of some models, hence the name.

Our estimator proceeds by obtaining an estimate, $\hat{\eta}$, of $\eta$. All
nuisance parameters are estimated using an ensemble regression approach known
as the Super Learner (Van der Laan et al., 2007; van der Laan & S. Dudoit &
A.W. van der Vaart, 2006), which creates an optimally weighted combination of
algorithms. Then, the pseudo-outcome $D_{\hat{\eta}}(O_{i})$ is computed for
all individuals in the sample, and a regression of the pseudo-outcome on
baseline covariates $V$ is fitted. Following the approach of Luedtke and van
der Laan (2016), we also use Super Learner in fitting the regression of
$D_{\hat{\eta}}(O)$ on $V$ and use the fitted values to obtain predictions
$\hat{\mathsf{B}}(v)$. The group of individuals at risk for a harmful indirect
effect is identified as those indices $i$ with
$\hat{\mathsf{d}}(v_{i})=\mathds{1}\\{\hat{\mathsf{B}}(v_{i})>0\\}$.

We use cross-fitting (Klaassen, 1987; Zheng and van der Laan, 2011;
Chernozhukov et al., 2016) throughout the estimation process in all fits. Let
${\cal V}_{1},\ldots,{\cal V}_{J}$ denote a random partition of data with
indices $i\in\\{1,\ldots,n\\}$ into $J$ prediction sets of approximately the
same size such that $\bigcup_{j=1}^{J}{\cal V}_{j}=\\{1,\ldots,n\\}$. For each
$j$, the training sample is given by ${\cal
T}_{j}=\\{1,\ldots,n\\}\setminus{\cal V}_{j}$. $\hat{\eta}_{j}$ denotes the
estimator of $\eta$, obtained by training the corresponding prediction
algorithm using only data in the sample ${\cal T}_{j}$, and $j(i)$ denotes the
index of the validation set which contains observation $i$. We then use these
fits, $\hat{\eta}_{j(i)}(O_{i})$ in computing $D_{\hat{\eta}_{j(i)}}(O_{i})$.
Likewise, regressions of the pseudo-outcome $D_{\hat{\eta}_{j(i)}}(O_{i})$ on
$V_{i}$ are trained within each training sample, and the final estimate is
computed by predicting in the corresponding validation data set.

##  4 Estimating the indirect effect under a hypothetical treatment decision
$\mathsf{d}(v)$

As stated in the Introduction, our goal is to identify a subset of the
population that is predicted to have a harmful indirect effect, and for whom
treatment may be carefully considered or reconsidered. Even though our goal is
not to develop a treatment rule, in some situations it may be important to
assess the population effects that would be observed if the function
$\mathsf{d}$ were used as a treatment rule. To that end, in this section we
estimate the hypothetical population interventional indirect effect if we were
to use $\mathsf{d}(V)$ to assign treatment to each individual. Specifically,
we define the total effect of implementing $\mathsf{d}(V)$ as
$\mathsf{E}(Y_{\mathsf{d}}-Y_{0})$, and decompose it in terms of direct and
indirect effects as

$\mathsf{E}(Y_{\mathsf{d}}-Y_{0})=\underbrace{\mathsf{E}(Y_{\mathsf{d},G_{\mathsf{d}}}-Y_{\mathsf{d},G_{0}})}_{\text{indirect
effect}}+\underbrace{\mathsf{E}(Y_{\mathsf{d},G_{0}}-Y_{0,G_{0}})}_{\text{direct
effect}}.$

The population interventional indirect effect can be identified under the
sequential randomization assumptions and positivity in Section 3.1 as:

$\mathsf{E}(Y_{\mathsf{d},G_{\mathsf{d}}}-Y_{\mathsf{d},G_{0}})=\\\
\int\mathsf{b}(\mathsf{d}(v),z,m,w)\mathsf{q}(z\mid\mathsf{d}(v),w)\\{\mathsf{p}(m\mid\mathsf{d}(v),w)-\mathsf{p}(m\mid
0,w)\\}\mathsf{p}(w)\mathrm{d}\nu(w,z,m).$ (3)

Because we do not have the true, unknown $\mathsf{d}$, but instead have an
estimate of it, $\hat{\mathsf{d}}$, we estimate the population indirect effect
$\mathsf{E}(Y_{\hat{\mathsf{d}},G_{\hat{\mathsf{d}}}}-Y_{\hat{\mathsf{d}},G_{0}})$
that would be observed if our estimate were implemented. We use the one-step
estimation approach described in Díaz et al. (2019), and which is based on
solving the EIF estimating equation. Specifically, we let our estimator be
defined as

$\hat{\theta}=\frac{1}{n}\sum_{i=1}^{n}\\{D_{\hat{\eta}}^{\hat{\mathsf{d}},\hat{\mathsf{d}}}(O_{i})-D_{\hat{\eta}}^{\hat{\mathsf{d}},0}(O_{i})\\},$

We use a cross-fitted version of this estimator, as described in Section 3.2.
The variance of this estimator can be estimated as the sample variance of the
EIF. Theorem 5 of van der Laan and Luedtke (2015) proves that plugging in
$\hat{\mathsf{d}}$ estimated from the data results in asymptotically linear
estimation of its effect, even though the same data were used to estimate
$\hat{\mathsf{d}}$ and to assess its effect.

The R code to implement this cross-fitted estimator is available:
blindedforreview.

##  5 Identifying those with predicted harmful indirect effects

We now apply our proposed approach to identify the subgroup of boys with a
predicted harmful indirect effect of voucher receipt on adolescent alcohol use
through aspects of the school, neighborhood, and social environments. In this
case, the average treatment effect of Section 8 housing voucher receipt on
long-term risk of adolescent alcohol use among boys is beneficial,
contributing to decreased risk, whereas the indirect effect through mediators
related to the school and neighborhood environments and instability of the
social environment is estimated to contribute to an increased risk of
subsequent adolescent alcohol use. There is evidence that improvements in
measures of the school and neighborhood environment (e.g., neighborhood
poverty, school rank, school poverty) may actually negatively impact a
subgroup of boys, possibly via increased risk of suspensions/expulsions and
less social support (Rudolph et al., 2020), which in turn, may increase risk
of a variety of negative mental health and substance use outcomes (Rudolph et
al., 2020).

If one is interested in optimizing the overall effect, one could apply
existing methods to identify the subgroup who may be harmed by the
intervention (Luedtke and van der Laan, 2016; Robins, 2004; Zhao et al.,
2012). However, instances where the overall effect is beneficial, but the
indirect effect is harmful—like this one—motivate our proposed approach to
identify a subgroup predicted to experience an unintended, harmful indirect
effect, even in the presence of a predicted beneficial overall effect. In this
application, we 1) identify the subgroup of boys who would be predicted to
have a harmful indirect effect and thus, for whom voucher receipt may not be
recommended, and 2) estimate the hypothetical indirect effect and hypothetical
total effect had the identified subgroup not received the intervention.

A limitation of using so-called “black-box” algorithms proposed in Section 3.2
for this task is that we are left without knowing which characteristics are
important in predicting whether an indirect effect will be positive or
negative. Consequently, we also use a recently proposed adaptive lasso
(Bahamyirou et al., 2020) to learn less-than-optimal, but interpretable
$\mathsf{d}$.

###  5.1 Data and Analysis

MTO was a large-scale experiment in which families living in low-income public
housing could sign-up to be randomized to receive a Section 8 housing voucher,
which is a form of housing assistance that subsidizes rent on the private
market, allowing families to move out of public housing. We restrict to
adolescent boys who were surveyed at the final follow-up timepoint in the
Boston, Chicago, New York City, and Los Angeles sites (N=2,100, rounded sample
size). We consider baseline covariates at the individual, family, and
neighborhood level ($W$, study site, age, race/ethnicity, number of family
members, previous problems in school, enrolled in special class for gifted and
talented students, parent is high school graduate, parent marital status,
parent work status, receipt of AFDC/TANF, whether any family member has a
disability, perceived neighborhood safety, neighborhood satisfaction,
neighborhood poverty level, reported reasons for participating in MTO,
previous number of moves, previous application for Section 8 voucher),
baseline randomized intervention status of received voucher or not ($A$,
binary 1/0), whether or not the family used the voucher to move ($Z$, binary
1/0), mediating variables representing aspects of the school and neighborhood
environments and instability of the social environment over the 10-15 years of
follow-up ($M$, school rank, student-to-teacher ratio, % students receiving
free or reduced-price lunch, % schools attended that were Title I, number of
schools attended, number of school changes within the year, whether or not the
most recent school was in the same district as the baseline school, number of
moves, neighborhood poverty; all weighted over the duration of follow-up), and
the long-term outcome of past-month alcohol use at the final timepoint, when
the children were adolescents ($Y$, binary 1/0). For the purposes of this
illustration, we use just one imputed dataset. As recommended for all MTO
analyses, we use individual-level weights that account for the sampling of
children within families, assignment ratios, and loss-to-follow-up
(Sanbonmatsu et al., 2011).

We first apply the estimation approach in Section 3.2 to identify the subgroup
$\\{i\in\\{1,\ldots,n\\}:\hat{\mathsf{d}}(V_{i})=0\\}$ of boys with a
predicted harmful indirect effect of voucher receipt on adolescent alcohol use
through aspects of the school, neighborhood, and social environments. We then
apply the estimation procedure in Section 4 to estimate indirect effects:
$\mathsf{E}(Y_{\hat{\mathsf{d}},G_{\hat{\mathsf{d}}}}-Y_{\hat{\mathsf{d}},G_{0}})$,
$\mathsf{E}(Y_{1,G_{1}}-Y_{1,G_{0}})$, and total effects:
$\mathsf{E}(Y_{\hat{\mathsf{d}},G_{\hat{\mathsf{d}}}}-Y_{0,G_{0}})$,
$\mathsf{E}(Y_{1,G_{1}}-Y_{0,G_{0}})$. We use five folds in cross-fitting, and
include: a simple mean model, main-terms generalized linear model, lasso,
multivariate adaptive regression splines, and extreme gradient boosted
machines in our Super Learner ensemble. Finally, we report the estimated
indirect effects and total effects using the alternative, interpretable
$\hat{\mathsf{d}}$ estimated using adaptive lasso.

###  5.2 Results

Figure 1 shows: the typical population interventional indirect effect (PIIE,
$\mathsf{E}(Y_{1,G_{1}}-Y_{1,G_{0}})$) estimates and total effect (PITE,
$\mathsf{E}(Y_{1,G_{1}}-Y_{0,G_{0}})$) estimates not using any
individualization, which are labeled as “no individualization”, and the PIIE
and PITE estimates using $\hat{\mathsf{d}}$ to hypothetically assign boys to
receive the voucher who were not predicted to have a harmful experience
through mediators of the school, neighborhood, and social environments
(learned via Super Learner), denoted “superlearner estimation”. Using Super
Learner to identify the subgroup predicted to experience a harmful indirect
effect and then not giving the voucher to that subgroup would be expected to
reduce the otherwise harmful indirect effect down to 0.0004 increased risk of
subsequent adolescent alcohol use (95% CI: -0.0690, 0.0698) as opposed to
0.0432 increased risk (95% CI: -0.1009, 0.1873). It would also be expected to
slightly decrease the beneficial total average effect of voucher assignment to
0.0224 reduced risk of subsequent adolescent alcohol use (-0.0224, 95% CI:
-0.1500, 0.1053) as opposed to 0.0683 reduced risk (95% CI: -0.2882, 0.1053).

Applying the adaptive lasso to estimate $\hat{d}$, we find that those whose
parent graduated from high school and whose race is not white, black, nor
Hispanic/Latino would be predicted to not have a harmful indirect effect. The
PIIE and PITE estimates using this $\hat{\mathsf{d}}$ are denoted “lasso
estimation”.

Figure 1: Population interventional indirect effects and total effects by rule
type. All results were approved for release by the U.S. Census Bureau,
authorization number CBDRB-FY21-ERD003-004.

##  6 Conclusions

We proposed a nonparametric approach to identify a subgroup who may have
unintended harmful indirect effects even in the presence of beneficial overall
effects. We then proposed a nonparametric estimator to estimate the
hypothetical population interventional indirect and total effects were one to
make the decision not to treat the subgroup predicted to have harmful indirect
effects. Our cross-fitted estimators solve the efficient influence function,
so are robust and can incorporate machine learning in model fitting.

This work was motivated by surprising results from a large, housing policy
intervention (the Moving to Opportunity Study) in which boys whose families
were randomized to receive a housing voucher had slightly reduced risk of
subsequent alcohol use as adolescents (total effect) but the indirect pathway
from voucher receipt to that positive effect operating through aspects of the
school and social environments was harmful. This harmful indirect effect could
reflect processes like boys whose families moved with the vouchers going to
schools with higher academic performance but where they were more likely to
face discriminatory practices, like suspension and expulsion, and more social
isolation (Rudolph et al., 2020). It is possible that some may wish to avoid
such negative processes, prompting the desire to identify those at risk for
them. In this example, we do so and find that not giving the intervention to
this subgroup would have indeed reduced the harmful predicted indirect effect
point estimate, though the confidence intervals are wide and overlapping.

In terms of identification and estimation, our proposal builds on both the
causal mediation and individualized optimal treatment regime literatures.
However, we caution against interpreting $\hat{\mathsf{d}}$ as a decision
rule, understanding that, at least in some instances, deciding whether or not
someone would benefit from an intervention/treatment is necessarily more
nuanced than the mediation mechanisms and total effects we consider. Future
work could focus on incorporating additional complexities and nuances into the
existing tools used in personalized medicine to reflect the multiple
mechanisms by which treatments or interventions may affect numerous, and
sometimes competing, outcomes among individuals, recognizing heterogeneities
(anticipated or not) at each stage in the process.

## References

  * Avin et al. (2005) Avin, C., I. Shpitser, and J. Pearl (2005).  Identifiability of path-specific effects.
  * Bahamyirou et al. (2020) Bahamyirou, A., M. E. Schnitzer, E. H. Kennedy, L. Blais, and Y. Yang (2020).  Doubly robust adaptive lasso for effect modifier discovery.  arXiv preprint arXiv:2011.12746.
  * Chernozhukov et al. (2016) Chernozhukov, V., D. Chetverikov, M. Demirer, E. Duflo, C. Hansen, et al. (2016).  Double machine learning for treatment and causal parameters.  arXiv preprint arXiv:1608.00060.
  * Clampet-Lundquist et al. (2011) Clampet-Lundquist, S., K. Edin, J. R. Kling, and G. J. Duncan (2011).  Moving teenagers out of high-risk neighborhoods: How girls fare better than boys.  American Journal of Sociology 116(4), 1154–89.
  * Díaz et al. (2019) Díaz, I., N. S. Hejazi, K. E. Rudolph, and M. J. van der Laan (2019).  Non-parametric efficient causal mediation with intermediate confounders.  arXiv preprint arXiv:1912.09936.
  * Díaz et al. (2018) Díaz, I., O. Savenkov, and K. Ballman (2018).  Targeted learning ensembles for optimal individualized treatment rules with time-to-event outcomes.  Biometrika 105(3), 723–738.
  * Didelez et al. (2006) Didelez, V., A. P. Dawid, and S. Geneletti (2006).  Direct and indirect effects of sequential treatments.  In Proceedings of the Twenty-Second Conference on Uncertainty in Artificial Intelligence, pp. 138–146. AUAI Press.
  * Kennedy (2020) Kennedy, E. H. (2020).  Optimal doubly robust estimation of heterogeneous causal effects.  arXiv preprint arXiv:2004.14497.
  * Kessler et al. (2014) Kessler, R. C., G. J. Duncan, L. A. Gennetian, L. F. Katz, J. R. Kling, N. A. Sampson, L. Sanbonmatsu, A. M. Zaslavsky, and J. Ludwig (2014).  Associations of housing mobility interventions for children in high-poverty neighborhoods with subsequent mental disorders during adolescence.  JAMA 311(9), 937–948.
  * Klaassen (1987) Klaassen, C. A. (1987).  Consistent estimation of the influence function of locally asymptotically linear estimators.  The Annals of Statistics, 1548–1562.
  * Kling et al. (2007) Kling, J. R., J. B. Liebman, and L. F. Katz (2007).  Experimental analysis of neighborhood effects.  Econometrica 75(1), 83–119.
  * Laber et al. (2014) Laber, E. B., K. A. Linn, and L. A. Stefanski (2014).  Interactive model building for q-learning.  Biometrika 101(4), 831–847.
  * Lok (2016) Lok, J. J. (2016).  Defining and estimating causal direct and indirect effects when setting the mediator to specific values is not feasible.  Statistics in medicine 35(22), 4008–4020.
  * Lok (2019) Lok, J. J. (2019).  Causal organic direct and indirect effects: closer to baron and kenny.  arXiv preprint arXiv:1903.04697.
  * Luedtke and van der Laan (2016) Luedtke, A. R. and M. J. van der Laan (2016).  Super-learning of an optimal dynamic treatment rule.  The international journal of biostatistics 12(1), 305–332.
  * MacKinnon et al. (2000) MacKinnon, D. P., J. L. Krull, and C. M. Lockwood (2000).  Equivalence of the mediation, confounding and suppression effect.  Prevention science 1(4), 173–181.
  * Moodie et al. (2012) Moodie, E. E., B. Chakraborty, and M. S. Kramer (2012).  Q-learning for estimating optimal dynamic treatment rules from observational data.  Canadian Journal of Statistics 40(4), 629–645.
  * Murphy (2003) Murphy, S. (2003).  Optimal dynamic treatment regimes.  Journal of the Royal Statistical Society: Series B 65(2), 331–66.
  * Nabi et al. (2018) Nabi, R., P. Kanki, and I. Shpitser (2018).  Estimation of personalized effects associated with causal pathways.  In Uncertainty in artificial intelligence: proceedings of the… conference. Conference on Uncertainty in Artificial Intelligence, Volume 2018\. NIH Public Access.
  * Pearl (2000) Pearl, J. (2000).  Causality: models, reasoning, and inference.  Cambridge University Press. ISBN 0 521(77362), 8.
  * Petersen et al. (2006) Petersen, M. L., S. E. Sinisi, and M. J. van der Laan (2006).  Estimation of direct causal effects.  Epidemiology, 276–284.
  * Qian and Murphy (2011) Qian, M. and S. A. Murphy (2011).  Performance guarantees for individualized treatment rules.  Annals of statistics 39(2), 1180.
  * Robins (2003) Robins, J. (2003).  Discussion of “optimal dynamic treatment regimes" by susan a. murphy.  Journal of the Royal Statistical Society: Series B 65(2), 355–66.
  * Robins (2004) Robins, J. M. (2004).  Optimal structural nested models for optimal sequential decisions.  In D. Y. Lin and P. J. Heagerty (Eds.), Proceedings of the second seattle Symposium in Biostatistics, pp. 189–326. Seattle, Washington: Springer Science+ Business Media.
  * Rubin and van der Laan (2007) Rubin, D. and M. J. van der Laan (2007).  A doubly robust censoring unbiased transformation.  The international journal of biostatistics 3(1).
  * Rudolph et al. (2020) Rudolph, K. E., C. Gimbrone, and I. Diaz (2020).  Helped into harm: Mediation of the effect of a housing voucher intervention on long-term mental health and substance use outcomes in adolescent boys.  arXiv preprint arXiv:2006.07708.
  * Rudolph et al. (2018) Rudolph, K. E., O. Sofrygin, N. M. Schmidt, R. Crowder, M. M. Glymour, J. Ahern, and T. L. Osypuk (2018).  Mediation of neighborhood effects on adolescent substance use by the school and peer environments.  Epidemiology 29(4), 590–598. PMCID: PMC5987191.
  * Rudolph et al. (2017) Rudolph, K. E., O. Sofrygin, W. Zheng, and M. J. Van Der Laan (2017).  Robust and flexible estimation of stochastic mediation effects: A proposed method and example in a randomized trial setting.  Epidemiologic Methods 7(1).
  * Sanbonmatsu et al. (2011) Sanbonmatsu, L., L. F. Katz, J. Ludwig, L. A. Gennetian, G. J. Duncan, R. C. Kessler, E. K. Adam, T. McDade, and S. T. Lindau (2011).  Moving to opportunity for fair housing demonstration program: Final impacts evaluation.
  * Schmidt et al. (2017) Schmidt, N. M., M. M. Glymour, and T. L. Osypuk (2017).  Adolescence is a sensitive period for housing mobility to influence risky behaviors: An experimental design.  Journal of Adolescent Health 60(4), 431–437.
  * Shpitser and Sherman (2018) Shpitser, I. and E. Sherman (2018).  Identification of personalized effects associated with causal pathways.  In Uncertainty in artificial intelligence: proceedings of the… conference. Conference on Uncertainty in Artificial Intelligence, Volume 2018\. NIH Public Access.
  * van der Laan and Luedtke (2015) van der Laan, M. J. and A. R. Luedtke (2015).  Targeted learning of the mean outcome under an optimal dynamic treatment rule.  Journal of causal inference 3(1), 61–95.
  * van der Laan and Petersen (2008) van der Laan, M. J. and M. L. Petersen (2008).  Direct effect models.  The international journal of biostatistics 4(1).
  * Van der Laan et al. (2007) Van der Laan, M. J., E. C. Polley, and A. E. Hubbard (2007).  Super learner.  Statistical applications in genetics and molecular biology 6(1).
  * van der Laan & S. Dudoit & A.W. van der Vaart (2006) van der Laan & S. Dudoit & A.W. van der Vaart, M. (2006).  The cross-validated adaptive epsilon-net estimator.  Statistics & Decisions 24(3), 373–395.
  * VanderWeele and Tchetgen Tchetgen (2017) VanderWeele, T. J. and E. J. Tchetgen Tchetgen (2017).  Mediation analysis with time varying exposures and mediators.  Journal of the Royal Statistical Society. Series B, Statistical Methodology 79(3), 917.
  * VanderWeele et al. (2014) VanderWeele, T. J., S. Vansteelandt, and J. M. Robins (2014).  Effect decomposition in the presence of an exposure-induced mediator-outcome confounder.  Epidemiology (Cambridge, Mass.) 25(2), 300.
  * Vansteelandt and Daniel (2017) Vansteelandt, S. and R. M. Daniel (2017).  Interventional effects for mediation analysis with multiple mediators.  Epidemiology (Cambridge, Mass.) 28(2), 258.
  * Zhang et al. (2012) Zhang, B., A. A. Tsiatis, M. Davidian, M. Zhang, and E. Laber (2012).  Estimating optimal treatment regimes from a classification perspective.  Stat 1(1), 103–114.
  * Zhao et al. (2012) Zhao, Y., D. Zeng, A. J. Rush, and M. R. Kosorok (2012).  Estimating individualized treatment rules using outcome weighted learning.  Journal of the American Statistical Association 107(499), 1106–1118.
  * Zhao et al. (2015) Zhao, Y.-Q., D. Zeng, E. B. Laber, and M. R. Kosorok (2015).  New statistical learning methods for estimating optimal dynamic treatment regimes.  Journal of the American Statistical Association 110(510), 583–598.
  * Zheng and van der Laan (2017) Zheng, W. and M. van der Laan (2017).  Longitudinal mediation analysis with time-varying mediators and exposures, with application to survival outcomes.  Journal of causal inference 5(2).
  * Zheng and van der Laan (2011) Zheng, W. and M. J. van der Laan (2011).  Cross-validated targeted minimum-loss-based estimation.  In Targeted Learning, pp. 459–474. Springer.
  * Zheng and van der Laan (2012) Zheng, W. and M. J. van der Laan (2012).  Targeted maximum likelihood estimation of natural direct effects.  The international journal of biostatistics 8(1), 1–40.


# Dynamics of horizontal buoyant jets.

V. P. Goncharov<EMAIL_ADDRESS>A. M. Obukhov Institute of Atmospheric
Physics RAS, 109017 Moscow, Russia

(August 27, 2024)

###### Abstract

A single-layer model is used to construct the theory of straight-line jets. In
addition to stationary states, the evolution of such jets admits transverse
pulsations. The theory predicts that the period for warm jets pulsations is
longer than the period of inertial oscillations caused by the Earth’s
rotation, while for cold jets pulsations, it is shorter. Thus, only warm jets
can have a noticeable effect on the synoptic range of the wind-speed spectrum,
generating into it additional spectral peaks. Within such an approach, every a
single jet is described by a self-similar compactly localized solution that
can be considered as a constant-shear band.

††preprint: AIP/123-QED

##  I Introduction

Jet streams in the Earth’s atmosphere are fast flowing, narrow air currents
produced by two factors: the atmospheric heating and the influence of the
Coriolis force caused by the planet rotation (rhi94, ; val06, ; hh13, ). Being
one of the elements in the atmospheric dynamics, jet streams can be
responsible for strong long-period fluctuations of the horizontal wind. In
particular, the presence of additional peaks in the synoptic range of the
spectrum of horizontal wind speed we can consider as a confirmation of this
effect.

The main spectral properties of the horizontal wind velocity established by
Van der Hoven (vdh57, ) are shown in Fig. 1. This spectrum has three well-
defined peaks. One of them called the synoptic maximum is at the period about
$4\,days$. The next peak corresponds to half a day, which is close to the
period of inertial oscillations caused by the Earth’s rotation. And the last
known as the microscale maximum locates at a period of $1\,min$. The frequency
range between the synoptic and microscale peaks contains almost no energy and
is known as the spectral gap.

Beginning with the classic work of Van der Hoven (vdh57, ) and up to the
present day, there has been accumulating enough experimental evidence that the
macro-meteorological spectral range may also include spectral peaks for
fluctuations with periods exceeding inertial. Examples of such observations
are presented in Kang and Won (kw15, ; kw16, )(see also the references quoted
therein).

A reasonable explanation of additional spectral peaks and the macro-
meteorological spectral maximum can be obtained within shallow water theory,
assuming that the jets are horizontal narrow flows of one fluid at the bottom
of another and evolve under the action of gravity and the Coriolis force
created by the rotation of the Earth.

Depending on the ratio between the fluid density inside such a jet and the
background density of fluid around it, its evolution has three alternative
scenarios. The joint action of the Coriolis and the Archimedean (buoyancy)
forces either causes transverse pulsations of this jet, or leads to its
narrowing and increasing in height, or makes it monotonically expand with
decreasing in altitude.

Figure 1: Power spectral density of horizontal wind speed by Van der Hoven
(vdh57, ).

It is noteworthy that similar behavior resulting from effects of the buoyancy
and rotation is detected for vortex structures considered within the active
layer model (del03, ; gp12, ; gp13a, ; gp13b, ; gp15, ; gp16, ). These
structures look like drops or fingers that are inside a homogeneous fluid and
have elliptic contact with a horizontal, rigid bottom. Being by solutions of a
conservative dynamical system, they can rotate and change their semi-axes, and
hence, the contact area and the height without changing the volume, total
angular momentum as well as other integrals of motion.

The drop- or finger-like structures are self-similar solutions and inherent to
many of the so-called active layer models, which describe a depth-averaged
flow of incompressible fluid with a horizontally-nonuniform density in the
lower (active) layer. Application of such models provides a way to understand
the different real processes initiated by the Rayleigh-Taylor instability. In
addition, since these structures produce space-time singularities responsible
for the power-law tails in the short-wave range of the spectrum, their study
provides a key to the understanding of strong turbulence (knnr07, ).

As practice has shown, using upward jets as structural elements of the theory,
these models can quite adequately describe the vertical mixing, transport
processes, and turbulence in astrophysical and geophysical flows (del03, ;
gp12, ; gp13a, ; gp13b, ; gp15, ; gp16, ). The results obtained in these works
inspire a hope that horizontal jet currents can also get a description within
the active layer models. One argument in favor of this is that, in the
limiting case, when one of the ellipse semi-axes in the base of drop tends to
zero, the vortex drop-like structures in these models can transform into
infinitely thin horizontal jets.

The generation of oscillations by horizontal jets due to their transverse
pulsations deserves a more detailed discussion and will be presented in the
next sections. Here, we note only that pulsational periods of warm jets should
be longer than the period of inertial oscillations and therefore lie in the
synoptic range. This fact can be argued with purely qualitative
considerations.

Let’s imagine that a fluid fragment has a density contrast with ambient fluid
and is under the influence of gravity and the Coriolis forces. If the rotation
is absent, the action of gravity leads to only two regimes. The fluid fragment
either sinks to the bottom or rises upwards, depending on whether or not it is
heavier or lighter than the ambient fluid. Because of irreversibility, both
effects can be treated as instability. And conversely, if we eliminate
gravity, then the action of the Coriolis force can cause only oscillations.
Thus, we can conclude that both the instabilities and oscillations are regimes
equally accessible for the active layer models.

According to the theory of dynamic systems, such symbiosis of regimes implies
a specific topology of dynamical trajectories in the phase space of the
system. In this space, stable regimes corresponding to oscillations evolve
along closed trajectories, while the unstable along the non-closed. In doing
so, the zones of stability and instability separate each other by closed loops
known as separatrices. It is due to the active layer models have separatrices
the oscillation periods caused by gravity and by the Earth’s rotation can
undergo significant changes in response to a small increase in the amplitude
of these oscillations.

Although the present paper limits to study only oscillatory regimes for
horizontal jets, the Rayleigh-Taylor instability can also make a notable
impact on the spectrum of wind speed fluctuations. The latter mechanism is of
importance when studying large-scale turbulence, sometimes called a structural
because its structural elements are narrow upward jets. In particular, a
similar approach was used to explain a scaling effect in the spectra of the
Sun’s supergranulation (gp14, ). In this work, the supergranulation was
considered as a collective effect produced by a statistical ensemble of upward
jets.

This paper is structured as follows. In section II, we discuss the active
layer model in the context of its following application to study horizontal
jet streams. In section III, this model is used to derive the master equations
for describing straight-line jets with transverse pulsations. In section IV,
we analyze self-similar solutions and give examples for jet profiles that
could provide these kinds of solutions. The pulsational dynamics, impacting
not only transverse parameters of jets but also components of their velocity,
are discussed in section V. One of the goals of this section is to show that
the period of inertial fluctuations resulting from the Coriolis effect is the
shortest period for pulsations of warm jets. Possible nonlinear mechanisms to
increase this period are also discussed here. In section VI, we summarize our
results. The appendix provides some details omitted in section V when
integrating the nonlinear differential equation.

##  II Active layer model

We consider a two-dimensional field model whose evolution is governed by the
equations

$\displaystyle\frac{du}{dt}-fv=-\frac{1}{2}h\tau_{x}-\tau h_{x},$ (1)
$\displaystyle\frac{dv}{dt}+fu=-\frac{1}{2}h\tau_{y}-\tau h_{y},$ (2)
$\displaystyle\frac{dh}{dt}+h\left(u_{x}+v_{y}\right)=0,\quad\frac{d\tau}{dt}=0.$
(3)

The notation is as follows: $x$, $y$ are the Cartesian horizontal coordinates,
$d/dt$ is the total temporal derivative, subscripts denote the partial
derivatives, and $u$, $v$ are the depth-averaged components of the horizontal
velocity.

Eqs. (1)–(3) describe the depth-averaged flow of incompressible fluid in an
boundary layer of thickness $h(x,y,t)$ in shallow-water approximation, on
assuming of hydrostatic balance and of rotation with the constant angular
velocity $\frac{1}{2}f$ about the vertical axis. The members $fv$ and $fu$ on
the left side of Eqs. (1), (2) implies the components of the Coriolis
acceleration. The physical specification of the variable $\tau$ depends on
what assumption about the density distribution is used in the model.

If we suppose that the boundary layer rests on an impermeable horizontal
bottom and is bounded above by a free surface (rip93, ), then the variable
$\tau$ has the meaning of the total buoyancy and is defined as

$\tau=\textsl{g}(1+\Delta\varrho/\varrho_{0}),$ (4)

where g is the gravity, and $\Delta\varrho$ is the density deviation from the
background value $\varrho_{0}$. According to this definition, the total
buoyancy $\tau$, along with the thickness $h$, is a strictly positive function
of horizontal coordinates and time. The formulation based on Eqs. (1)–(3) with
such definition of buoyancy is known as the thermal rotating shallow-water
model (TRSW) (see Zeitlin (ztl18, ), sec. 14 and references therein).

Note that in the case when $\Delta\varrho=0$, so that $\tau=\textsl{g}$ is
constant, Eqs. (1)–(3) reduce to the usual shallow water equations describing
motion under the influence of gravity g and the Coriolis force.

In the atmosphere where the free boundary condition is hardly applicable, it
is relevant to refuse it by transforming the TRSW model into the so-called
active layer model (del03, ; gp12, ; gp13a, ; gp13b, ; gp15, ). For this, all
we need to do is add one more layer of an incompressible fluid with the
density $\varrho_{0}$ so that it fills all the half-space above the bottom
layer up to infinity.

If both of these fluids are under the action of gravity and Coriolis forces,
and the horizontally-nonuniform density jump $\Delta\varrho$ between them is
small, so that $\Delta\varrho/\varrho_{0}\ll 1$, then in such models (del03, ;
gp12, ; gp13a, ; gp13b, ; gp15, ) the shallow-water approximation, the depth-
averaging, and the requirement of hydrostatic balance leads us again to Eqs.
(1)–(3). But now, instead of (4), these models use the relation

$\tau=\textsl{g}\Delta\varrho/\varrho_{0},$ (5)

which gives any sign for $\tau$. Thus, in the active layer model, the variable
$\tau$ has the meaning of the relative buoyancy, rather than merely buoyancy
as in the condition (4), where it takes only positive values.

When density variations are induced only by temperature ones $\Delta T$ and
are linearly connected, the relative buoyancy can be computed as

$\tau=-\textsl{g}\Delta T/T_{K},$ (6)

where $T_{K}$ is the background absolute (Kelvin) temperature. This
parametrization allows taking into account both heating and cooling effects
(see as examples (and84, ; rip99, )).

As with any non-dissipative system, the active-layer model has symmetries
which imply the existence of conservation laws. Among them are the total
energy

$H=\frac{1}{2}\int h\left(u^{2}+v^{2}+h\tau\right)dxdy,$ (7)

the total angular momentum $S=M+fI/2$, and the quantity

$G=\frac{1}{4}\left(\frac{dI}{dt}\right)^{2}+\left(M+2f^{-1}H\right)^{2}$ (8)

that remains unchanged due to the scaling invariance (gp14, ). Here the
angular momentum $M$ and the moment of inertia $I$ are described by the
expressions:

$M=\int h(xv-yu)dxdy,\quad I=\int h\left(x^{2}+y^{2}\right)dxdy.$ (9)

Additionally, the model leaves invariant the integral

$C=\int h\left(F_{1}(\tau)+qF_{2}(\tau)\right)dxdy,$ (10)

where $F_{1,2}$ are arbitrary functions of the buoyancy $\tau$, and $q$ is the
potential vorticity, which is defined as

$q=h^{-1}\left(f+v_{x}-u_{y}\right),$ (11)

and obeys the equation

$\frac{dq}{dt}=\frac{1}{2h}\left(h_{x}\tau_{y}-h_{y}\tau_{x}\right).$ (12)

##  III Straight-line jets

Next, the active layer model will be applied to study so-called straight-line
jets on the assumption that their dynamics includes only the changes in
direction and transverse pulsations as shown in Fig. 2. To formulate the
motion equations for these jets in a more easy-to-use form, we introduce,
instead of $x$ and $y$, new coordinates

$s=x\cos\theta+y\sin\theta,\quad n=y\cos\theta-x\sin\theta.$ (13)

Here $s$ and $n$ are the longitudinal and transverse coordinates accompanying
the jet’s angular motion, and $\theta(t)$ is the angle of its slope relative
to $x$-axis.

Figure 2: Scheme of a straight-line jet flow, which is oriented along the axis
$s$ and pulsates transversely to it.

In the coordinates $s$, $n$, new velocities are defined as

$U=\frac{ds}{dt},\quad V=\frac{dn}{dt},$ (14)

and are connected with the old ones by the relations

$\displaystyle
u=U\cos\theta-V\sin\theta-\left(s\sin\theta+n\cos\theta\right)\theta_{t},$
(15) $\displaystyle
v=U\sin\theta+V\cos\theta+\left(s\cos\theta-n\sin\theta\right)\theta_{t}.$
(16)

Taking into account these formulas, after some simplifications, we find that
in the new coordinates the equations of motion (1)–(3) take the form

$\displaystyle\frac{dU}{dt}-2V\theta_{t}-s\theta_{t}^{2}-n\theta_{tt}-f\left(V+s\theta_{t}\right)=-h_{s}\tau-\frac{1}{2}h\tau_{s},$
(17)
$\displaystyle\frac{dV}{dt}+2U\theta_{t}-n\theta_{t}^{2}+s\theta_{tt}+f\left(U-n\theta_{t}\right)=-h_{n}\tau-\frac{1}{2}h\tau_{n},$
(18) $\displaystyle
h_{t}+\partial_{s}\left(Uh\right)+\partial_{n}\left(Vh\right)=0,$ (19)
$\displaystyle\tau_{t}+U\tau_{s}+V\tau_{n}=0.$ (20)

We are looking for solutions that can be physically interpreted as homogeneous
(along the $s$-axis) jets whose dynamics and therefore the variables $h$,
$\tau$, $U$, $V$ depend only on time and the transverse coordinate $n$. As the
analysis shows, for such jet-currents, the angular velocity $\theta_{t}$ must
either be absent $(\theta_{t}=0)$ or constant $(\theta_{t}=-f)$. For both
these cases the equations of motion (17)–(20) can be rearranged to the form

$\displaystyle\tau_{t}+V\tau_{n}=0,$ (21) $\displaystyle U_{t}+VU_{n}-\omega
V=0,$ (22) $\displaystyle h_{t}+\partial_{n}\left(Vh\right)=0,$ (23)
$\displaystyle V_{t}+VV_{n}+\omega U=-h_{n}\tau-\frac{1}{2}h\tau_{n},$ (24)

where the parameter $\omega$ takes two values: $\omega=f$ if $\theta_{t}=0$,
and $\omega=-f$ if $\theta_{t}=-f$.

It is easy to verify that, as before, Eqs. (21)–(24) conserve the integrals
(7) and (10), but only in a one-dimensional form

$\displaystyle H=\frac{1}{2}\int h\left(U^{2}+V^{2}+h\tau\right)dn,$ (25)
$\displaystyle C=\int
h\bigl{(}F_{1}\left(\tau\right)+qF_{1}\left(\tau\right)\bigr{)}dn,$ (26)

where now the role of the potential vorticity $q$ is played by the quantity

$q=h^{-1}\left(\omega-U_{n}\right),$ (27)

which according to (22), (23) obeys the equation

$\frac{dq}{dt}=q_{t}+Vq_{n}=0.$ (28)

As for the rest of the laws of conservation, they turn out to be unavailable
in the one-dimensional version because reducing the dimensionality leads to
the loss of certain symmetries.

##  IV Self-similar solutions

By limiting to the class of self-similar solutions, when examining Eqs.
(21)–(24) we will assume that the buoyancy $\tau$ depends on $n$ and $t$ in a
self-similar manner as

$\tau=\hat{\tau}(z),\quad z=n\beta^{-1},$ (29)

thus the new independent variable $z$ is the self-similarity variable, and
$\beta=\beta(t)$ is a function of time.

The assumption of self-similarity allows us to make a few simplifications.
Indeed, by using the ansatz (29) in Eq. (21), we at once obtain the expression
for the transverse velocity

$V=-\tau_{t}/\tau_{n}=\beta_{t}z.$ (30)

Substitution of this result into Eqs. (22), (23) gives us three more
constraints

$h=\beta^{-1}\hat{h}(z),\quad U=\omega\alpha z,\quad\alpha_{t}=\beta_{t},$
(31)

where $\hat{h}(z)$ is some function which along with $\hat{\tau}(z)$ is to be
found. In besides, the condition imposed on the longitudinal velocity $U$
allows us to interpret the obtained solution as a shear band, which from
outside looks like a flow composing of two counter jets.

Finally, utilizing the obtained results (30), (31) and applying the method of
separation of variables, from Eq. (24) it follows that, on the one hand,

$\beta_{tt}+f^{2}\alpha-\kappa\beta^{-2}=0,$ (32)

and on the other hand, the profile functions $\hat{h}(z)$ and $\hat{\tau}(z)$
obey the equation

$\partial_{z}\bigl{(}\hat{h}^{2}\hat{\tau}\bigr{)}+2\kappa z\hat{h}=0.$ (33)

Here $\kappa$ is an arbitrary constant known as the separation constant.

Since the velocities $V$ and $U$ depend linearly on z, to be physically
reasonable, the solutions of Eq. (33) must be entirely localized on a finite
interval of $z$ and vanishing outside of it. Besides, as the quantity $h$ is
strictly positive, due to the invariance with respect to the replacement
$z\rightarrow-z$, Eq. (33) can have only even solutions. Therefore, without
loss of generality, we may seek its solutions on the interval $-1\leq z\leq 1$
in assuming they are symmetric about $z=0$ and subject to the boundary
conditions

$\hat{h}|_{z=\pm 1}=0,\quad\hat{\tau}|_{z=\pm 1}=0.$ (34)

The physical meaning of the separation constant $\kappa$ becomes clear if
paying attention to the fact that it connects with two integrals of the
problem

$\mathfrak{U}=\frac{1}{2}\int h^{2}\tau dn,\quad I=\int n^{2}hdn.$ (35)

One of them $\mathfrak{U}$ is the potential energy, and the other $I$ is known
as the total momentum inertia. In the self-similar case, both integrals can be
presented as

$\displaystyle\mathfrak{U}=\beta^{-1}\mathfrak{U}_{0},\quad I=\beta^{2}I_{0},$
(36)
$\displaystyle\mathfrak{U}_{0}=\frac{1}{2}\int_{-1}^{1}\hat{h}^{2}\hat{\tau}dz,\quad
I_{0}=\int_{-1}^{1}z^{2}\hat{h}dz,$ (37)

where $\mathfrak{U}_{0}$ and $I_{0}$ are constants functionally depended on
profile functions $\hat{h}$ and $\hat{\tau}$. Multiplying Eq. (33) on $z$ and
integrating over the interval $-1\leq z\leq 1$, one can show that

$\kappa=\mathfrak{U}_{0}/I_{0}.$ (38)

Another property can be directly derived from Eq. (33) if rewriting it as

$\hat{\tau}=\frac{2\kappa}{\hat{h}^{2}}\int_{z}^{1}z\hat{h}dz,$ (39)

Thus, as follows from this relation, the signs of the separation constant
$\kappa$ and the relative buoyancy $\hat{\tau}$ always coincide.

To construct solutions for the profile functions $\hat{h}$ and $\hat{\tau}$,
we will consider a generating function $\psi(z)$, which has the following
properties. First, being defined in the interval $-1\leq z\leq 1$, it changes
the sign of its derivative $d\psi/dz$ in the point $z=0$, and second, near the
points $z=\pm 1$ it behaves as

$\psi\approx\left(1-z^{2}\right)^{1+\gamma},$ (40)

where the parameter $\gamma$ satisfies the inequality $0\leq\gamma\leq 1$.

Then for any such function $\psi(z)$, solutions of Eq. (33) can be presented
in the following form

$\hat{h}=-\frac{1}{z}\frac{d\psi}{dz},\quad\hat{\tau}=2\kappa\frac{\psi}{\hat{h}^{2}}.$
(41)

It is easy to check that under the condition (40), the solutions (41) possess
all the necessary properties. They provide positivity of the function
$\hat{h}$ and turn it together with $\hat{\tau}$ into zero at points $z=\pm 1$
according to the power laws

$\hat{h}\propto\left(1-z^{2}\right)^{\gamma},\quad\hat{\tau}\propto\left(1-z^{2}\right)^{1-\gamma}.$
(42)

Despite the functional arbitrariness in the choice of the generating function
$\psi$, all jet-streams with the same parameter $\kappa$, have the identical
transverse dynamics. This fact follows immediately from Eq. (32), since under
fixed $\omega$ the evolution of the self-similarity variable $\beta(t)$
depends only on $\kappa$.

Two pairs of the solutions with different profiles, but with the same
parameter $\kappa$, are presented in Figs. 3, 4.

Figure 3: The normalized profiles in cross-section of the jet, which are
generated by the function (44). The solid line shows the thickness $\hat{h}/a$
and the dashed line displays the buoyancy $\hat{\tau}a/\kappa$.

The first of them includes the profiles

$\hat{h}=a\left(1-z^{2}\right)^{1/2},\quad\hat{\tau}=\frac{2\kappa}{3a}\left(1-z^{2}\right)^{1/2},$
(43)

corresponding to the generating function

$\psi=\frac{a}{3}\left(1-z^{2}\right)^{3/2},$ (44)

where $a>0$ is some arbitrary constant.

Figure 4: The normalized profiles in cross-section of the jet, which are
generated by the function (47). The solid line shows the thickness $\hat{h}/a$
and the dashed line displays the buoyancy $\hat{\tau}a/\kappa$.

The other pair consists of the profiles

$\displaystyle\hat{h}=ae^{-z^{2}}\sqrt{1-z^{2}},$ (45)
$\displaystyle\hat{\tau}=\frac{\kappa}{a}\frac{e^{z^{2}}}{\sqrt{1-z^{2}}}\left(1-\frac{D\left(\sqrt{1-z^{2}}\right)}{\sqrt{1-z^{2}}}\right),$
(46)

generated by the function

$\psi=\frac{a}{2}e^{-z^{2}}\left(\sqrt{1-z^{2}}-D\left(\sqrt{1-z^{2}}\right)\right).$
(47)

where $D$ is the Dawson function.

##  V Pulsations of jets

Within the framework of this model, the pulsational dynamics of jet structures
involves the changes not only in their transverse parameters but also in the
longitudinal component of velocity. Completely presented by Eqs. (30)–(32),
the jet dynamics is described by two variable $\alpha$ and $\beta$ which
according to (30) are related as

$\alpha=\beta-c,$ (48)

where $c$ is an integration constant.

The elimination of $\alpha$ allows us to formulate the motion equation (32) in
the Hamiltonian form

$\beta_{t}=H_{\lambda},\quad\lambda_{t}=-H_{\beta}$ (49)

with the Hamiltonian $H$, defined by the expression

$H=\frac{1}{2}\left(\lambda^{2}+f^{2}\beta^{2}\right)+\frac{\kappa}{\beta}-cf^{2}\beta.$
(50)

The constant $c$ can be fixed by assuming that the oscillatory system has a
stable stationary state in which it obeys the conditions

$H_{\lambda}=0,\quad H_{\beta}=0,\quad H_{\beta\beta}>0.$ (51)

Whence, supposing that this state is achieved at a value of $\beta=\beta_{0}$,
we find the two constraints

$c=\beta_{0}-\frac{\kappa}{f^{2}\beta_{0}^{2}},\quad\frac{2\kappa}{f^{2}\beta_{0}^{3}}>-1.$
(52)

In order to reduce the number of parameters, we make the problem
dimensionless, by using $\beta_{0}$ and $f^{-1}$ as the characteristic scales
of length and time. Then after some algebra, from (50) we get

$\displaystyle\eta_{t}^{2}=2\left(\mathfrak{H}-\mathfrak{U}\right),$ (53)
$\displaystyle\mathfrak{H}=\frac{1}{2}-k+\frac{H}{f^{2}\beta_{0}^{2}},\quad\mathfrak{U}=\frac{1}{2}\left(\eta-1\right)^{2}\left(1+\frac{k}{\eta}\right),$
(54)

Here the variable $\eta$ and the parameter $k$ are dimensionless quantities
given as

$\eta=\frac{\beta}{\beta_{0}},\quad k=\frac{2\kappa}{f^{2}\beta_{0}^{3}},$
(55)

Eq. (53) can be considered as the equation of motion of a point particle in
the potential $\mathfrak{U}(\eta)$. As the analysis shows (see Fig. 5), the
model has two oscillatory regimes, depending on the parameter $k$, whose sign
coincides with the sign of the temperature deviation $\Delta T$.

Figure 5: The domains of existence for oscillatory regimes. Their boundaries
are lines of equal periods $T=0$, $T=2\pi$, $T=\infty$. The dark gray region
corresponds to warm jets with $-1<k<0$, and the light gray one – cold jets
with $k>0$.

An oscillatory hierarchy in the theory of straight-line jets is arranged in
such a way that only warm jets with $-1<k<0$ can provide long-period
pulsations with $T>2\pi$. Although cold jets with $k>0$ also supports
pulsations, their periods can only be less than $2\pi$. Moreover, as shown in
Fig. 6, in contrast to cold jets, the phase trajectories of oscillations for
warm jets are surrounded by separatrices.

Figure 6: A sketch of the potential and phase portrait for the model in the
case $\kappa<0$. The solid line marks the separatrix loop, and the oscillatory
domain within it is colored gray.

This fact has direct relevance to the subject-matter of this paper, since any
separatrix loop implies the existence of the finite amplitudes in the system,
at approaching to which the oscillation period can become arbitrarily large.
For the model under consideration, in particular, this provides a mechanism
through which the oscillation period due to the Earth’s daily rotation can
undergo significant changes under the influence of nonlinearity.

Assuming that the quantity $\eta$ oscillates in the range
$\eta_{2}\leq\eta\leq\eta_{3}$, Eq. (53) can be transformed into

$\eta_{t}^{2}=\eta^{-1}\left(\eta-\eta_{1}\right)\left(\eta-\eta_{2}\right)\left(\eta_{3}-\eta\right),$
(56)

where, in force of the condition $-1<k<0$, the right side roots obey to
inequalities:

$0<\eta_{1}<\eta_{2}<1<\eta_{3}.\\\ $ (57)

On integration of Eq. (56) (see the Appendix for details), using the
parameters $\eta_{2}$ and $\eta_{3}$ instead of $\mathfrak{H}$ and $k$ is
preferable. In this case, it is necessary to have in mind that

$\displaystyle\mathfrak{H}=\frac{1}{2}\frac{(1-\eta_{2})^{2}(1-\eta_{3})^{2}}{1-\eta_{2}\eta_{3}},$
(58) $\displaystyle
k=\frac{\eta_{2}\eta_{3}(\eta_{2}+\eta_{3}-2)}{1-\eta_{2}\eta_{3}},\quad\eta_{1}=\frac{2-\eta_{2}-\eta_{3}}{1-\eta_{2}\eta_{3}}.$
(59)

As shown in the Appendix, by introducing the notations

$p=1-\frac{\eta_{2}}{\eta_{3}},\quad m=\frac{\eta_{1}}{\eta_{2}-\eta_{1}}.$
(60)

the solution of (56) can be presented in parametric form

$\displaystyle\eta=\frac{\eta_{2}}{1-p\sin^{2}\varphi},$ (61) $\displaystyle
t=2\sqrt{(1-p)(1+m)}\;\Pi(p;\varphi|-pm),$ (62)

where $\varphi$ is the phase, and $\Pi$ is the incomplete elliptic integral of
the third kind.

Since the solution (61), (62) describes oscillations, their period $T$ can be
found from the condition $T=t|_{\varphi=\pi/2}$, which leads us to the
expression

$T=4\sqrt{(1-p)(1+m)}\;\Pi(p;-pm).$ (63)

It should be underlined this quantity is bounded from below $T\geq 2\pi$.
That, in turn, means that if a nonlinear system with $-1<k<0$ is under the
action of Coriolis force, inertial oscillations $T=2\pi$ should manifest
themselves in it as a minimum possible period. According to (60), this is
exactly the value that is reached when $p\rightarrow 0$ or $m\rightarrow 0$.

As may be concluded from the conditions (57), oscillatory regimes for warm
jets are realized with the accordance of the inequalities

$\frac{2}{1+\eta_{2}}\leq\eta_{3}\leq 2-\eta_{2}.$ (64)

In the space of parameters $\eta_{2}$ and $\eta_{3}$, they define a narrow
existence domain (see Fig. 5), at whose boundaries the oscillation period $T$
and the parameter $k$ take the following values

$\displaystyle T\Bigr{|}_{\eta_{3}=2-\eta_{2}}=2\pi,\quad
k\Bigr{|}_{\eta_{3}=2-\eta_{2}}=0,$ (65) $\displaystyle
T\Bigr{|}_{\eta_{3}=2/(1+\eta_{2})}=\infty,\quad
k\Bigr{|}_{\eta_{3}=2/(1+\eta_{2})}=-\frac{2\eta_{2}^{2}}{1+\eta_{2}}.$ (66)

For comparison, note that for cold jets at the upper boundary of the domain,
where $\eta_{3}=1/\eta_{2}$, these quantities have the following values

$T\Bigr{|}_{\eta_{3}=1/\eta_{2}}=0,\quad
k\Bigr{|}_{\eta_{3}=1/\eta_{2}}=\infty.$ (67)

From the last Eqs. (64)–(66) we can drawn a number of conclusions. Firstly,
the upper boundary in Fig. 5 corresponds to the simplest case of linear
(harmonic) oscillations, which develop against the background of neutral
thermal stratification $k=0$. Thus, in this regime, the equality

$a_{+}=a_{-},\quad a_{+}=\eta_{3}-1,\quad a_{-}=1-\eta_{2}.$ (68)

for amplitude deviations $a_{+}$ and $a_{-}$ from the stationary state is
supported exactly.

In the whole rest of the domain, the violation of this equality is quite weak.
For each fixed $\eta_{2}$, the most asymmetry in oscillatory amplitudes is
achieved at the lower boundary, where parameter $k$ takes the threshold value
$k=-2\eta_{2}^{2}/(1+\eta_{2})$. In the point $\eta_{2}=\sqrt{2}-1$, the
amplitude difference $\Delta=a_{-}-a_{+}$ reaches a peak value equal to
$\Delta\approx 0{.}17$.

The only quantity that is substantially changing in the domain is the period
$T$. As shown in Fig. 7, this effect manifests itself as the exponential
growth of $T$.

Figure 7: The period $T$ at $\eta_{2}=1-\sqrt{2}$ as a function of the
deviation of the amplitude parameter $\eta_{3}$ from the threshold value
$\sqrt{2}$.

It occurs at once the amplitudes $\eta_{2}$ and $\eta_{3}$ become sufficiently
close to thresholds. These values are reached on the lower boundary of the
domain of existence and related one another as

$\eta_{3}=\frac{2}{1+\eta_{2}}.$ (69)

Such behavior indicates that the phase trajectory of the oscillation is near
the separatrix. As a result, the period of such oscillations can be very
different from the inertial period $T=2\pi$, which is the minimum permissible
for the given system.

One more scenario leading to disturbances, whose temporal durations are very
large in comparison with the inertial period, occurs when their phase
trajectories turn to be exactly on the separatrix. In this case the parameters
$\eta_{1}$ and $\eta_{2}$ coincide, $\eta_{3}$ satisfies the equality (69),
and Eq. (56) can be transformed as

$\eta_{t}^{2}=\eta^{-1}\left(\eta-\eta_{2}\right)^{2}\left(\frac{2}{1+\eta_{2}}-\eta\right).$
(70)

Such solutions are called separatrix ones. They initially increase, reaching
at some time the maximum value $\eta=2/(1+\eta_{2})$, and then decrease to an
asymptotic level $\eta=\eta_{2}>0$.

Based on this scenario, to integrate Eq. (70) explicitly, we assume that, in
the course of evolution, the system moves along the separatrix loop. In doing
so let it leave the saddle point at $t=-\infty$, stop at $t=0$, and return to
the saddle point at $t=+\infty$. In this case, as a result, we get

$\displaystyle
t=T_{s}\ln\left(\frac{\sqrt{\zeta}-T_{s}\sqrt{1-\zeta}}{\sqrt{\zeta}+T_{s}\sqrt{1-\zeta}}\right)-2\arccos\sqrt{\zeta},$
(71) $\displaystyle
T_{s}=\sqrt{\frac{\eta_{2}(1+\eta_{2})}{(1-\eta_{2})(2+\eta_{2})}},\quad\eta=\frac{2\zeta}{1+\eta_{2}}.$
(72)

The solution describes the temporal evolution during which the jet first
expands and then narrows, restoring its original width at
$t\rightarrow\infty$. It should be noted that the separatrix solutions have
only one key parameter $k$. As shown in Fig. 8, all other parameters depend on
it, including $T_{s}$, which plays the role of a characteristic time.

Figure 8: The time scale $T_{s}/(2\pi)$ and amplitude parameters $\eta_{2}$ an
$\eta_{3}$ as a functions of $k$ for separatrix solutions.

##  VI Discussion and conclusions

This paper deals with studying the dynamics of horizontal jet-streams and
addresses at least two issues. The first one concerns the principal
possibility of using the active layer model to achieve this goal. And the
second relates to the applicability of the obtained results to some problems
of geophysical fluid dynamics. The explanation of additional peaks in the
macrometeorological range of the horizontal wind speed spectrum is one such
problem.

The straight-line jet theory developed within the active layer model, in
addition to stationary states, assumes two more dynamical possibilities for
jets: transverse pulsations and rotation at a constant angular velocity. The
theory also predicts that warm jets pulsate with periods longer than the
period of inertial oscillations caused by the Coriolis force, while cold jets
pulsate with periods shorter than inertial oscillations. So, only warm jets
which satisfy the condition

$-1<k<0,$ (73)

where $k$ defined in (55) is a dimensionless factor of the thermal
stratification, can contribute to the macrometeorological spectral range.

According to the results of the theory, both jet velocity components
(longitudinal one $U$ and transverse $V$) are expressed through the variable
$\eta$ as follows

$V=\omega\eta_{t}\frac{n}{\eta},\quad
U=\omega\frac{n}{\eta}\left(\eta-1+k\right).$ (74)

The extremal values of these quantities are achieved at jet stream boundaries
$n=\pm\beta_{0}$ and have opposite signs at them.

In the stationary state when $\eta=1$, these expressions yield the formulas

$V=0,\quad U=\omega kn.$ (75)

which describe a uniform shear flow in a band of constant width $2\beta_{0}$.
Because the profile of the longitudinal velocity $U$ is asymmetrical relative
to $n=0$, such flow, in some sense, can be considered as a result of
superposing two opposite-directed jets.

Depending on the parameter $\omega=\pm f$, the band with shear flow either
remains immovable or rotates at the constant angular velocity $-f$. In the
first case, when $\omega=f$ and hence $k<0$, the flow in the immovable band
has positive vorticity

$\zeta=\omega-U_{n}=\omega(1-k).$ (76)

In the second case, when $\omega=-f$ and the band is rotating, its flow
vorticity becomes negative.

In order, relying on the jet theory, to deeper reach an insight into the
investigated phenomena and correctly interpret the observable data, it is
necessary to have an appropriate criterion. A suitable criterion can be
derived directly from Eq. (75) and the inequality (73) if the maximal
longitudinal velocity along with the half-width $\beta_{0}$ are considered as
typical values. Then this criterion will look like

$|k|=\frac{U_{0}}{f\beta_{0}}<1.$ (77)

Here $U_{0}$ is the modulus of the maximum value of the longitudinal velocity
reached at jet boundaries $n=\pm\beta_{0}$.

The criterion can be modified so that it will be able to take into account the
influence of both geometric and temperature effects on jet pulsations. To this
end, we need to express the stratification factor $k$ through the rest
parameters of the jet. Using Eqs. (43), (6) and applying them to axial values,
we get the relationships for the jet height $h_{0}$ and the buoyancy
$\tau_{0}$:

$h_{0}\tau_{0}=\frac{2\kappa}{3\beta_{0}},\quad\tau_{0}=-\textsl{g}\frac{\Delta
T_{0}}{T_{K}}.$ (78)

Here $\Delta T_{0}$ is the axial value of the temperature deviation from the
background absolute temperature $T_{K}$.

Using (78) and (55), the stratification factor $k$ can be presented as

$k=-3Bu,\quad Bu=\textsl{g}\frac{\Delta
T_{0}}{T_{K}}\frac{h_{0}}{f^{2}\beta_{0}^{2}},$ (79)

where $Bu$ is the so-called Burger number.

For instance, the conditions:

$\displaystyle g=9{.}8\,m\cdot s^{-2},\quad\beta_{0}=2\cdot 10^{5}\,m,\quad
h_{0}=10^{3}\,m,$ (80) $\displaystyle f=1{.}2\cdot 10^{-4}\,s^{-1},\quad\Delta
T_{0}/T_{K}=10^{-2},$ (81)

which are typical of the Earth’s atmosphere, in accordance with (79), lead us
to the equality $k=-0{.}51$. It means that in a stationary state, the
longitudinal velocity of the jet takes the value $U_{0}=12{.}25\,m\cdot
s^{-1}$.

If the amplitude of such a jet is pulsating in the range of $\eta_{2}$ and
$\eta_{3}$, these values must obey the constraint (59). In particular, if we
take $\eta_{3}=1{.}2$, then at $k=-0{.}51$ we find $\eta_{2}=0{.}7$. The jet
pulsations with these parameters are shown in Fig. 9. Their period calculated
in accordance with (63) amounts to $T=1{.}93\,days$.

Figure 9: Synoptic jet pulsations. Two-day pulsations correspond to a solid
line and seven-day pulsations to a dashed line.

As for the transverse velocity, it reaches its maximum at $\eta=1$, which
allows us to reach the expression

$V_{0}=U_{0}\frac{(1-\eta_{2})(\eta_{3}-1)}{\sqrt{1-\eta_{2}\eta_{3}}}.$ (82)

The calculations by this formula for the jet considered above give
$V_{0}=1{.}84\,m\cdot s^{-1}$.

For comparison, the dashed line in Fig. 9 presents pulsations for one more
jet. However, despite it has very close in amplitude parameters
$\eta_{2}=0{.}6666667$ and $\eta_{3}=1{.}2$, its period $T=6{.}9\,days$
differs noticeably. It means that near the separatrix, there is a high
sensitivity of the oscillation period to the jet parameters.

Obviously, in real geophysical systems, it is almost impossible that the
background parameters are stationary with such high accuracy. Perhaps that is
why it is quite problematic to experimentally observe regular oscillations
with periods of more than 2-3 days.

## Appendix A

Let’s consider Eq. (56) and enter the notations

$p=1-\frac{\eta_{2}}{\eta_{3}},\quad m=\frac{\eta_{1}}{\eta_{2}-\eta_{1}}.$
(83)

The solution of Eq. (56) will be searching in the parametric form

$\eta=\frac{\eta_{2}}{1-p\sin^{2}\varphi},\quad t=t\left(\varphi\right),$ (84)

where the phase $\varphi$ plays the role of the parameter.

Then for variable $t(\varphi)$ satisfies the equation

$t_{\varphi}^{2}=\frac{4(1-p)(1+m)}{\left(1-p\sin^{2}\varphi\right)^{2}\left(1+pm\sin^{2}\varphi\right)},$
(85)

which has the solution

$t=2\sqrt{(1-p)(1+m)}\Pi(p;\varphi|-pm),$ (86)

where $\Pi$ is the incomplete elliptic integral of the third kind.

Since the positive parameters $\eta_{1}$, $\eta_{2}$, and $\eta_{3}$ are
arranged as $\eta_{1}<\eta_{2}\leq\eta_{3}$, the parameters $p$ and $m$
satisfy to inequalities

$0\leq p<1,\quad m\geq 0.$ (87)

###### Acknowledgements.

This work was equally supported by the Russian Foundation for Basic Research
(Project No. 18-05-00831), by the Presidium of the Russian Academy of Sciences
(Program of Nonlinear Dynamics in Mathematical and Physical Sciences), and by
the Russian Science Foundation (Project No. 19-17-00248).

## References

  * (1) P. B. Rhines, “Jets,” Chaos 4, 313–339 (1994).
  * (2) G. K. Vallis, Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-scale Circulation, (Cambridge University Press, 2006).
  * (3) J. Holton and G. Hakim, An Introduction to Dynamic Meteorology, Academic Press (Elsevier Science, 2013).
  * (4) I. Van der Hoven, “Power spectrum of horizontal wind speed in the frequency range from 0.0007 to 9000 cycles per hour,” Journal of Meteorology 14, 160–164 (1957).
  * (5) S. Kang and H. Won, “Intra-farm wind speed variability observed by nacelle anemometers in a large inland wind farm,” J. Wind Eng. Ind. Aerodyn. 147, 164–175 (2015).
  * (6) S. Kang and H. Won, “Spectral structure of 5 year time series of horizontal wind speed at the Boulder Atmospheric Observatory,” J. Geophys. Res. Atmos. 121, 946–967 (2016).
  * (7) P. J. Dellar, “Common Hamiltonian structure of the shallow water equations with horizontal temperature gradients and magnetic fields,” Physics of Fluids 15, 292–297 (2003).
  * (8) V. P. Goncharov and V. I. Pavlov, “Blow-up instability in shallow water flows with horizontally-non-uniform density,” JETP Lett. 96, 427–431 (2012).
  * (9) V. P. Goncharov and V. I. Pavlov, “Structural elements of collapses in shallow water flows with horizontally nonuniform density,” JETP 117, 754–763 (2013).
  * (10) V. P. Goncharov and V. I. Pavlov, “Simple model of the Rayleigh-Taylor instability, collapse, and structural elements,” Phys. Rev. E 88, 023002 (2013).
  * (11) V. P. Goncharov and V. I. Pavlov, “Algebraic instability in shallow water flows with horizontally nonuniform density,” Phys. Rev. E 91, 043004 (2015).
  * (12) V. P. Goncharov and V. I. Pavlov, “Pulsating jet-like structures in magnetized plasma,” Physics of Plasmas 23, 082117 (2016).
  * (13) E. A. Kuznetsov, V. Naulin, A. H. Nielsen, and J. J. Rasmussen, “Effects of sharp vorticity gradients in two-dimensional hydrodynamic turbulence,” Physics of Fluids 19, 105110 (2007).
  * (14) V. P. Goncharov and V. I. Pavlov, “Whether the Sun’s supergranulation possesses a scaling?” JETP Lett. 99, 317 (2014).
  * (15) P. Ripa, “Conservation laws for primitive equations models with inhomogeneous layers,” Geophys. Astrophys. Fluid Dynam. 70, 85–111 (1993).
  * (16) V. Zeitlin, Geophysical Fluid Dynamics: Understanding (almost) Everything with Rotating Shallow Water Models, (Oxford University Press, 2018).
  * (17) D. T. T. Anderson, “On equatorial dynamics, mixed layer physics and sea-surface temperature,” Tellus, Ser. A 36, 278–291 (1984).
  * (18) P. Ripa, “On the validity of layered models of ocean dynamics and thermodynamics with reduced vertical resolution,” Dynam. Atmos. Oceans 29, 1–40 (1999).