# Time-dependent treatment of cosmic-ray spectral steepening due to turbulence
driving

Martin Pohl University of Potsdam, Institute of Physics and Astronomy, D-14476
Potsdam, Germany  DESY, D-15738 Zeuthen, Germany

###### Abstract

Cosmic-ray acceleration at non-relativistic shocks relies on scattering by
turbulence that the cosmic rays drive upstream of the shock. We explore the
rate of energy transfer from cosmic rays to non-resonant Bell modes and the
spectral softening it implies. Accounting for the finite time available for
turbulence driving at supernova-remnant shocks yields a smaller spectral
impact than found earlier with steady-state considerations. Generally, for
diffusion scaling with the Bohm rate by a factor $\eta$, the change in
spectral index is at most $\eta$ divided by the Alfvénic Mach number of the
thermal sub-shock. For $M_{\mathrm{A}}\lesssim 50$ it is well below this
limit. Only for very fast shocks and very efficient cosmic-ray acceleration
the change in spectral index may reach $0.1$. For standard SNR parameters it
is negligible. Independent confirmation is derived by considering the
synchrotron energy losses of electrons: if intense nonthermal multi-keV
emission is produced, the energy loss, and hence the spectral steepening, is
very small for hadronic cosmic rays that produce TeV-band gamma-ray emission.

Gamma-ray astronomy — cosmic rays — supernova remnants — spectral index

††journal: ApJ

##  1 Introduction

Bell et al. (2019) studied cosmic-ray acceleration at nonrelativistic shocks
using a tensor expansion of the Vlasov equation, and they found a steepening
of the cosmic-ray spectrum at supernova remnants (SNR) arising from energy
transfer from the cosmic rays to turbulence in the precursor of the shocks.
Assuming the turbulence in question is nonresonant, the so-called Bell mode
(Bell, 2004), and taking an estimate for its magnetic-field energy density at
the saturation level,

$U_{\delta B}\approx\frac{v_{\mathrm{sh}}}{2c}\,U_{\mathrm{cr}},$ (1)

where $v_{\mathrm{sh}}$ denotes the shock speed and $U_{\mathrm{cr}}$ the
energy density in cosmic rays immediately upstream of the shock, they find a
spectral steepening by

$\Delta s|_{\mathrm{Bell}}\approx\frac{4}{\epsilon}\frac{U_{\delta
B}}{U_{\mathrm{cr}}}\approx\frac{2}{\epsilon}\frac{v_{\mathrm{sh}}}{c}.$ (2)

For young SNR this spectral change should be quite sizable. Here $\epsilon$ is
the fraction of turbulent energy density that is in the form of magnetic-field
fluctuations, $U_{\delta B}=\epsilon\,U_{\mathrm{wave}}$. For an incompressive
linear wave with phase speed $v_{\phi}<v_{\mathrm{A}}$, like Bell’s mode, we
have $\epsilon\gtrsim 0.5$. Clearly, in the nonlinear phase quite a bit of
energy is transferred to velocity fluctuations (Stroman et al., 2009) and
local displacements of electrons and ions (Kobzar et al., 2017), suggesting a
moderate variation in $\epsilon$. The modification of the bulk flow speeds
(Luo & Melrose, 2009; Kobzar et al., 2017; Weidl et al., 2019) is a nonlinear
transfer of wave energy to the plasma, as are secondary instabilities that
thrive on Bell modes (Bykov et al., 2011). These processes do not change
$\epsilon$, and so $\epsilon=0.5$ should be approximately valid for Bell modes
in the nonlinear stage.

The result conforms with intuition. The streaming of cosmic rays is the energy
supply of the turbulence. In the steady state the production rates of
turbulence and cosmic rays in the upstream region and at the shock must
balance the loss incurred by advection to the far-downstream region. The
advection speed is the same for the two. The energy density of the turbulence
is increased by compression at the shock, and so the fraction of the cosmic-
ray production power that is funneled to turbulence is a few times the ratio
of their upstream energy density, $U_{\delta B}/\epsilon U_{\mathrm{CR}}$. The
cosmic-ray spectrum reflects the balance between the energy gain and the
retention probability for each acceleration cycle (Bell, 1978). Energy
transfer to turbulence reduces the energy gain of cosmic rays, hence the
spectral steepening. Equation 2 then corresponds to the highest possible
effect, given the constraint $kr_{\mathrm{L}}\gg 1$ for Bell’s instability.

What of the uncertainties in the estimate? Does turbulence driving perhaps
impose termination of acceleration as opposed to a continuous spectral
steepening? Which fraction of the observed post-shock magnetic energy density
is actually carried by Bell’s mode, and what part is provided by other
processes operating at the shock (e.g. Giacalone & Jokipii, 2007)? Do the
waves at the shock fully sample the cosmic-ray energy loss or have the waves
lost energy already in the upstream region? And is the analytically estimated
saturation level really relevant? Is there enough time for Bell’s mode to grow
to very high amplitude in the precursor of a shock, before the plasma arrives
at the shock and the driving of waves stops? For simplicity we consider cosmic
rays with Lorentz factor $\Gamma_{\mathrm{cr}}$. The peak growth rate of the
mode scales with the proton gyrofrequency, $\Omega_{\mathrm{p}}$, as

$\gamma_{\mathrm{max}}\simeq\Omega_{\mathrm{p}}\,\frac{v_{\mathrm{sh}}\,N_{\mathrm{cr}}}{2\,v_{\mathrm{A}}\,N_{\mathrm{p}}}=\omega_{\mathrm{p,p}}\,\frac{v_{\mathrm{sh}}\,N_{\mathrm{cr}}}{2\,c\,N_{\mathrm{p}}},$
(3)

where $v_{\mathrm{A}}$ is the Alfvén speed, $\omega_{\mathrm{p,p}}$ stands for
the proton plasma frequency, and $N_{\mathrm{cr}}$ and $N_{\mathrm{p}}$ denote
the number density of cosmic rays and of ambient protons, respectively.

The cosmic-ray precursor has an extent that is determined by the ratio of the
diffusion coefficient, here written as multiple of the Bohm diffusion
coefficient, $\kappa=\eta r_{\mathrm{L}}c/3$, and the shock speed,
$v_{\mathrm{sh}}$. Strong driving of turbulence or cosmic-ray feedback can
moderately modify, and in situations with efficient particle acceleration
reduce, the extent. Magnetic-field amplification reduces the Larmor radius,
and Reville & Bell (e.g. 2013) find $\kappa$ to be well described by either
$\eta$ being a few for $r_{\mathrm{L}}$ as measured immediately upstream of
the shock or $\eta\lesssim 1$ with the far-upstream value of $r_{\mathrm{L}}$.
In this paper we shall use the former description and denote with $\delta B$
the amplitude of turbulent magnetic field very close to the shock, unless it
is explicitely written as $\delta B(r)$. Following the simulation results of
Reville & Bell (2013) we shall use $\eta=4$ as a fiducial value and consider
deviations from that in the discussion section. The shock-capture time needed
for the shock to cross a region as wide as the shock precursor is

$\tau_{\mathrm{sc}}\approx\frac{\kappa}{v_{\mathrm{sh}}^{2}}=\frac{\eta\,\Gamma_{\mathrm{cr}}\,c^{2}}{3\,\Omega_{\mathrm{p}}\,v_{\mathrm{sh}}^{2}}.$
(4)

The number of exponential growth cycles available for Bell’s mode then is
(Niemiec et al., 2008)

$N_{\mathrm{exp}}=\gamma_{\mathrm{max}}\tau_{\mathrm{sc}}\approx\frac{\eta\Gamma_{\mathrm{cr}}c^{2}N_{\mathrm{cr}}}{6\,v_{\mathrm{A}}v_{\mathrm{sh}}N_{\mathrm{p}}}=\frac{\eta
M_{\mathrm{A}}}{12}\frac{U_{\mathrm{cr}}}{U_{\mathrm{bulk}}},$ (5)

where $M_{\mathrm{A}}$ is the Alfvénic Mach number of the shock, and $U$
denotes the energy density in cosmic rays and in the bulk plasma flow
($=1/2\,\rho v_{\mathrm{sh}}^{2}$). Efficient cosmic-ray acceleration may
funnel $10\%$ of the bulk-flow energy into cosmic rays. Given our expectation
$\eta\approx 4$, it becomes clear that one needs Alfvénic Mach numbers of
$M_{\mathrm{A}}\approx 300$ or more to allow 10 exponential growth cycles of
the mode. Note that Bell’s mode is perceived to grow out of the fluctuation
spectrum in the interstellar medium or in the wind bubble of a core-collapse
supernova, and it is supposed to reach an amplitude much larger than that of
the large-scale field, $\delta B\gg B_{0}$, which realistically takes an
amplitude enhancement by much more than a factor thousand (or more than seven
exponential cycles). We note that even the notoriously noisy particle-in-cell
simulations require ten growth cycles, if not more, to reach significant
magnetic-field amplification (Riquelme & Spitkovsky, 2009; Stroman et al.,
2009; Gargaté et al., 2010; Kobzar et al., 2017). For a low upstream gas
density, $N_{\mathrm{p}}\lesssim 0.5\ \mathrm{cm}^{-3}$, and standard values
of the interstellar magnetic field, $B=7\ \mathrm{\mu G}$ (Jansson & Farrar,
2012), this is difficult to establish, on account of $v_{\mathrm{A}}\gtrsim
20\ \mathrm{km\,s^{-1}}$. In the environment of a core-collapse supernova the
magnetic field is less well known (Romanova & Owocki, 2015), but should
eventually scale inversely with the radius, $B\propto 1/r$ (Voelk & Forman,
1982), implying a spatially constant Alfvén speed in the free-wind zone.
Whatever the value of the Alfvén speed far upstream of the shock, significant
magnetic-field amplification in the cosmic-ray precursor will reduce it. For
high-density environments damping may strongly limit the growth of the mode
(Reville et al., 2007). Reaching a sufficient number of growth times for
magnetic-field amplification is certainly possible under specific conditions,
but not something one should expect for any SNR. If the conditions are met,
then for only a short period in the evolution of the remnant, because the
outer shock is likely too slow already at the end of the free-expansion phase.

Brose et al. (2020) demonstrated for Alfvénic turbulence that slow build-up of
turbulence in the upstream region will rapidly reduce the maximum energy, to
which an SNR can accelerate. In any case, if turbulence driving would
significantly soften the cosmic-ray spectrum, $U_{\mathrm{cr}}$ would fall off
with increasing energy. At very high energy only very few growth cycles would
be available and/or $\eta$ would be very large, both of which would stymie
further acceleration.

In this paper we present estimates of the energy loss of cosmic rays to
turbulence that account for the limited time available for the development of
the nonresonant streaming instability. In our discussion we shall use
parameter values that may represent the forward shocks of young SNRs, but in
principle the analysis applies to any fast non-relativistic shock at which
diffusive shock acceleration operates.

##  2 Time-dependent modeling of turbulence build-up

We shall explore two independent ways to estimate the energy transfer from
cosmic rays to turbulence. One will be a local consideration of turbulence
driving in the precursor, the other an integral assessment of the energy
transfer in the entire precursor. For simplicity we shall assume that all ions
in the cosmic rays and the background plasma are protons, i.e. hydrogen
nuclei. The non-resonant and broadband character of Bell’s mode permits
considering cosmic rays of a specific energy as a proxy for particles in a
wide energy band.

###  2.1 Local estimate

The growth rate of wave energy density is proportional to the energy density
in the unstable wave band,

$\dot{U}\simeq\int dk\
\gamma(k)\,\frac{B_{k}^{2}}{4\pi\epsilon}\lesssim\gamma_{\mathrm{max}}\,\frac{(\delta
B)^{2}}{4\pi\epsilon}=\gamma_{\mathrm{max}}\,\frac{2\,U_{\delta
B}}{\epsilon},$ (6)

where again $\epsilon$ is the fraction of the wave energy that is carried by
magnetic field, $\delta B$ is the magnetic wave amplitude, and $U_{\delta B}$
the magnetic-field energy density of the waves. Equation 6 indicates that most
of the energy transfer arises when the wave amplitude is high, meaning near or
at saturation. The expression is written as an upper limit that is
conservative for two reasons. First, it applies the peak growth rate to the
entire wave spectrum, including spectral bands in which only cosmic rays of
vastly different energy can drive Bell’s mode. Second, in the highly nonlinear
phase the growth rate, $\gamma_{\mathrm{max}}$, might be lower than it is for
moderate $\delta B$. Both effects can only strengthen the validity of our
upper limit for the growth rate of wave energy density. Inserting Eq. 3 gives

$\dot{U}\lesssim\frac{\omega_{\mathrm{p,p}}}{\epsilon}\frac{N_{\mathrm{cr}}}{N_{\mathrm{p}}}\frac{v_{\mathrm{sh}}}{c}U_{\delta
B}.$ (7)

Here $\omega_{\mathrm{p,p}}$ is the ion (proton) plasma frequency. To be noted
from eq. 7 is that the energy-density transfer rate is independent of the
amplitude of the large-scale magnetic field, $B_{0}$. It is evident that most
of the energy transfer occurs deep in the cosmic-ray precursor and close to
the shock, where both $N_{\mathrm{cr}}$ and $U_{\delta B}$ reach their
upstream peak values.

The energy gain by the waves imposes energy loss on the cosmic rays. We can
estimate the energy-loss time per cosmic-ray particle as

$\tau_{\mathrm{loss}}\simeq\frac{U_{\mathrm{cr}}}{\dot{U}}\gtrsim\frac{2\epsilon\Gamma_{\mathrm{cr}}}{\omega_{\mathrm{p,p}}}\frac{U_{\mathrm{bulk}}}{U_{\delta
B}}\frac{c^{3}}{v_{\mathrm{sh}}^{3}},$ (8)

which evidently is independent of the number density of cosmic rays. The only
property of the cosmic rays that enters Eq. 8 is their Lorentz factor.

We shall now compare the energy-loss time (Eq. 8) with the acceleration time
assuming diffusive shock acceleration. As the nonresonant streaming
instability operates only upstream of the shock, our comparison can ignore the
part of the acceleration time that is spent downstream. Particles spend
approximately half of their time in the upstream region, and so the effective
energy-loss time is twice that given in eq. 8,

$\tau_{\mathrm{loss,eff}}\simeq
2\tau_{\mathrm{loss}}\gtrsim\frac{4\epsilon\Gamma_{\mathrm{cr}}}{\omega_{\mathrm{p,p}}}\frac{U_{\mathrm{bulk}}}{U_{\delta
B}}\frac{c^{3}}{v_{\mathrm{sh}}^{3}}.$ (9)

Again expressing the diffusion coefficient of relativistic cosmic rays in Bohm
units, $\kappa=\eta c\,r_{\mathrm{L}}/3$, we rewrite equation (32) of Drury
(1991) as

$\tau_{\mathrm{acc}}=\frac{8\,\kappa}{v_{\mathrm{sh}}^{2}}=\frac{8\eta\Gamma_{\mathrm{cr}}}{3\Omega_{\mathrm{p}}}\frac{c^{2}}{v_{\mathrm{sh}}^{2}}.$
(10)

Following Bell (1978), the integral cosmic-ray spectrum reflects the balance
between acceleration and escape by advection to the far downstream region,
hence for the spectral index $s-1=\tau_{\mathrm{acc}}/\tau_{\mathrm{esc}}$.
Energy losses by driving turbulence increase the effective acceleration time
scale, leading to a softened spectrum

$\frac{dN(>E)}{dE}\approx\frac{\Delta N(>E)}{\Delta
E}\simeq\frac{1-s}{1-\frac{\tau_{\mathrm{acc}}}{\tau_{\mathrm{loss,eff}}}}\,\frac{N(>E)}{E}\
.$ (11)

For test particles at a strong shock in a hydrogen plasma one expects $s=2$.
The ratio of timescales in the denominator of eq. 11 must be less than unity,
otherwise acceleration is impossible. For a small ratio of acceleration time
and loss time the resulting change of spectral index is

$\displaystyle\Delta s$
$\displaystyle=\frac{s-1}{\frac{\tau_{\mathrm{loss,eff}}}{\tau_{\mathrm{acc}}}-1}\simeq(s-1)\frac{\tau_{\mathrm{acc}}}{\tau_{\mathrm{loss,eff}}}\quad\Rightarrow$
$\displaystyle\Delta s$
$\displaystyle\lesssim\frac{2\,(s-1)\,\eta}{3\,\epsilon}\frac{U_{\delta
B}}{U_{\mathrm{bulk}}}\frac{v_{\mathrm{sh}}}{c}\frac{\omega_{\mathrm{p,p}}}{\Omega_{\mathrm{p}}}.$
(12)

Here the proton gyrofrequency, $\Omega_{\mathrm{p}}$, derives from the scaling
of the acceleration time with the Bohm diffusion coefficient. Reville & Bell
(2013) noted that for efficient magnetic-field amplification the diffusion
coefficient is well described by the amplified field strength and $\eta\approx
4$. One can rewrite expression 12 using the Alfvénic Mach number,
$M_{\mathrm{A}}=v_{\mathrm{sh}}\omega_{\mathrm{p,p}}/c\Omega_{\mathrm{p}}$,
that is calculated with the amplitude of the amplified field without regard of
direction,

$\Delta
s\lesssim\frac{2\,(s-1)\,\eta\,M_{\mathrm{A}}}{3\,\epsilon}\frac{U_{\delta
B}}{U_{\mathrm{bulk}}}.$ (13)

###  2.2 Global assessment

We shall now conduct a global assessment of the cosmic-ray energy loss
incurred in the entire upstream region. The spatial profile of cosmic-ray
density in the so-called precursor is dominated by the homogeneous solution to
the spatial part of the transport equation,

$N_{\mathrm{cr}}\simeq
N_{\mathrm{cr,sh}}\exp\left(-\int_{r_{\mathrm{sh}}}^{r}\ \frac{dr^{\prime}\
v(r^{\prime})}{\kappa(r^{\prime})}\right),$ (14)

where $v(r)$ is the upstream flow speed measured in the shock rest frame. This
steady-state profile represents a balance between the advective flux toward
the shock and the diffusive flux away from it. The diffusive flux,
$-\kappa\,\partial N_{\mathrm{cr}}/\partial r$, also determines the cosmic-ray
current and hence the driving rate of Bell’s mode. If the cosmic-ray gradient
were weaker than its steady-state value, cf. Eq. 14, then advection would
dominate the spatial transport, and the density profile would evolve toward
the steady state; likewise for a steeper density profile. Equation 14 hence
represents a stable equilibrium and can be used in our subsequent
calculations.

The cosmic-ray density profile also determines the acceleration timescale,
because that depends on the average separation of upstream cosmic rays from
the shock, $\langle r-r_{\mathrm{sh}}\rangle$. In each acceleration cycle,
cosmic rays enjoy a relative energy gain on the order of $v_{\mathrm{sh}}/c$
and need a few times $\langle r-r_{\mathrm{sh}}\rangle/c$ for it. Hence the
acceleration time is a few times $\langle
r-r_{\mathrm{sh}}\rangle/v_{\mathrm{sh}}$. If $v$ and $\kappa$ were constant,
then

$\langle
r-r_{\mathrm{sh}}\rangle=\frac{\kappa}{v}=\frac{\kappa}{v_{\mathrm{sh}}}\,,$
(15)

and we recover the formula
$\tau_{\mathrm{acc}}\propto\kappa/v_{\mathrm{sh}}^{2}$ (Drury, 1991). If
$v/\kappa$ would significantly decrease at $r-r_{\mathrm{sh}}\gg\langle
r-r_{\mathrm{sh}}\rangle$, then we would arrive at approximately the same
conclusion, except that there might be some (weak) escape toward the far
upstream. If $v/\kappa$ would significantly decrease already close to the
shock, then $\langle r-r_{\mathrm{sh}}\rangle$ would be very large and may in
fact become unbound, in which case many of the freshly accelerated cosmic rays
escape to the far upstream, and the cosmic-ray spectrum would be very steep
(Brose et al., 2020). The rapid increase of volume with increasing $r$ would
contribute to the cut-off in the spectrum of confined particles (Ohira et al.,
2010).

In the steady state the cosmic-ray current would be unaffected and could still
drive the non-resonant mode (Reville et al., 2009), thus reducing $\kappa(r)$
and leading to a seemingly universal profile of the cosmic-ray precursor (Bell
et al., 2013). In essence, a significant increase of $\kappa$ with distance
from the shock would give ample time for turbulence driving, but would make
particle acceleration slow and inefficient, and hence impose a cut-off in the
spectrum. If that were the case, young SNRs could not have produced cosmic
rays of very high energy. The TeV-band detection of many SNRs suggests
otherwise, and so we will proceed with the assumption that the variation in
$\kappa(r)$ does not preclude a rapid decline of the cosmic-ray density in the
precursor. In other words, we consider cosmic-ray energies well below the
maximum the shock can presently provide. Studies of cosmic-ray-modified shocks
(Amato & Blasi, 2006) suggest that this approximation is good.

We shall now explicitly consider variations in the flow speed and hence write
eq. 3 with $v(r)$ instead of $v_{\mathrm{sh}}$. Using eq. 6 we integrate the
energy transfer rate over the entire precursor, which for a plane-parallel
shock means integration over $r$,

$\dot{E}_{\mathrm{tot}}\lesssim\frac{N_{\mathrm{cr,sh}}}{8\pi\,\epsilon\,c}\int_{r_{\mathrm{sh}}}^{\infty}dr\
\frac{\omega_{\mathrm{pp}}\,v(r)\,(\delta
B(r))^{2}}{N_{\mathrm{p}}\,\exp\left(\int_{r_{\mathrm{sh}}}^{r}dr^{\prime}\
\frac{v(r^{\prime})}{\kappa(r^{\prime})}\right)}.$ (16)

Continuity mandates that
$\omega_{\mathrm{pp}}/N_{\mathrm{p}}\propto\sqrt{v(r)}$. We can scale the
density and the flow speed to their values at the thermal sub-shock,
$N_{\mathrm{p,sh}}$ and $v_{\mathrm{sh}}$, and pull these quantities out of
the integral,

$\dot{E}_{\mathrm{tot}}\lesssim\frac{\omega_{\mathrm{pp,sh}}\,N_{\mathrm{cr,sh}}}{8\pi\,\epsilon\,c\,N_{\mathrm{p,sh}}\sqrt{v_{\mathrm{sh}}}}\int_{r_{\mathrm{sh}}}^{\infty}dr\
\frac{v(r)^{3/2}\,(\delta
B(r))^{2}}{\exp\left(\int_{r_{\mathrm{sh}}}^{r}dr^{\prime}\
\frac{v(r^{\prime})}{\kappa(r^{\prime})}\right)}.$ (17)

Note that $v(r)$ is measured in the sub-shock rest frame, and so non-linear
cosmic-ray feedback would lead to a positive gradient in $v(r)$ and a negative
gradient in $(\delta B(r))^{2}$ on account of compression, that partially
compensate each other.

The energy density in the magnetic turbulence likely increases toward the
shock on account of turbulence driving, and so $(\delta B(r))^{2}$ falls off
with increasing $r$. The cosmic-ray scattering rate is linked to the intensity
of the turbulence, hence $\kappa(r)$ will rise with increasing $r$. We shall
write the spatial profiles as

$(\delta B(r))^{2}=(\delta
B_{\mathrm{sh}})^{2}\,b(r)\qquad\kappa(r)=\kappa_{\mathrm{sh}}\,k(r).$ (18)

Note that $\delta B(r)$ may include small-$k$ turbulence that is driven by
cosmic rays of energies higher than that of the particles whose energy losses
we calculate. The diffusion coefficient increases with energy, and so does the
precursor length, implying that at any given location in the precursor a
particle sees turbulence that has been driven by cosmic rays of higher energy
further out in the precursor. Part of that turbulence has a scale commensurate
with or larger than the Larmor radius of the particle in question, which
provides additional scattering and possibly magnetic-field amplification, thus
reducing $\eta$ in our description, but it does not necessarily impose
additional energy loss on our particle.

The variation of $b(r)$ and $k(r)$ may be huge, as we observe evidence of
strong magnetic-field amplification. We noted above that $\kappa(r)$ should
not increase too quickly, otherwise the acceleration is inefficient and the
cosmic-ray spectrum cuts off.

We shall now use the argument of the exponential in eq. 17 as new variable of
integration,

$x=\int_{r_{\mathrm{sh}}}^{r}dr^{\prime}\,\frac{v(r^{\prime})}{\kappa(r^{\prime})},$
(19)

which leads to

$\displaystyle\dot{E}_{\mathrm{tot}}\lesssim\,$
$\displaystyle\frac{\omega_{\mathrm{pp,sh}}\,N_{\mathrm{cr,sh}}\,(\delta
B_{\mathrm{sh}})^{2}\,\kappa_{\mathrm{sh}}}{8\pi\,\epsilon\,c\,N_{\mathrm{p,sh}}}$
$\displaystyle\times\int_{0}^{\infty}dx\
\sqrt{\frac{v(x)}{v_{\mathrm{sh}}}}\,b(x)\,k(x)\,\exp(-x),$ (20)

where $b(x)$ and $k(x)$ are defined in eq. 18. In the quasilinear limit
$\kappa\propto 1/(\delta B)^{2}$ is approximately true, meaning
$b(x)k(x)\approx\mathrm{const}$, whereas for Bohmian scaling one expects
$b(x)k(x)\propto\sqrt{b(x)}$. In the general case, $b(x)k(x)$ is at most
constant, but more likely a declining function, because $b(x)$ is expected to
decline, and so the integral in eq. 20 yields a numerical factor close to
unity. Explicitly considering the spatial variation of $\delta B$ and of the
diffusion coefficient, $\kappa(r)$, thus leads to the same conclusion as the
phenomenological analysis described in sec. 2.1: a reasonable spatial
variation of the relevant parameters does not invalidate the calculations
presented in this section. The total energy-loss rate of cosmic rays in the
precursor then is

$\dot{E}_{\mathrm{tot}}\lesssim U_{\delta
B}\frac{N_{\mathrm{cr,sh}}}{N_{\mathrm{p}}}\frac{\eta\Gamma_{\mathrm{cr}}}{3\epsilon}\frac{c^{2}}{v_{\mathrm{A}}},$
(21)

where we again used Bohm scaling for the diffusion coefficient, $\kappa=\eta
r_{\mathrm{L}}c/3$. The escape flux to the far-downstream region,

$\dot{N}_{\mathrm{esc}}=-N_{\mathrm{cr,sh}}\frac{v_{\mathrm{sh}}}{4},$ (22)

is in the steady state compensated by the energy gain on account of
acceleration at the shock, and the level of balance determines the spectral
index, s,

$\dot{E}_{\mathrm{acc}}=\frac{m_{p}c^{2}\Gamma_{\mathrm{cr}}}{s-1}|\dot{N}_{\mathrm{esc}}|\
.$ (23)

The ratio of the two rates in eqs. 21 and 23 is a measure of the spectral
steepening,

$\Delta
s=(s-1)\frac{1}{\frac{\dot{E}_{\mathrm{acc}}}{\dot{E}_{\mathrm{tot}}}-1}\lesssim\frac{2\,(s-1)^{2}\,\eta\,M_{\mathrm{A}}}{3\,\epsilon}\frac{U_{\delta
B}}{U_{\mathrm{bulk}}}.$ (24)

As for eq. 12, we expanded this formula for weak losses to derive the last
expression. Apart from a factor $s-1$, we find the same level of spectral
steepening as in eq. 12. There we had used a formula for the acceleration time
scale that was strictly derived for test particles at a strong shock, i.e.
$s=2$, for which there is no difference between the two estimates for the
spectral steepening. We remind the reader that we integrated out spatial
variations in the cosmic-ray precursor and that all variables are supposed to
be measured immediately upstream of the thermal sub-shock.

##  3 Summary and Discussion

We calculated the energy-transfer rate from cosmic rays to non-resonant plasma
waves in the precursor of the forward shock of an SNR. Two different ways of
calculation led to essentially the same result for the softening of the
particle spectra that is imposed by that energy transfer. It can be related to
the energy-density ratio of amplified turbulent magnetic field and bulk plasma
flow. That ratio is at most the inverse square of the Alfvénic Mach number of
the thermal sub-shock, which we write with the full turbulent field amplitude,

$M_{\mathrm{A}}=\sqrt{\frac{U_{\mathrm{bulk}}}{U_{\delta
B}+U_{B_{0}}}}{\color[rgb]{0,0,0}=\sqrt{\frac{U_{\mathrm{bulk}}}{U_{\delta
B}}}\frac{1}{\sqrt{1+\frac{U_{B_{0}}}{U_{\delta B}}}}}.$ (25)

For a turbulently amplified magnetic field, $U_{\delta B}\gg U_{B_{0}}$, the
spectral steepening can hence be written either with the Alfvénic Mach number
of the sub shock or with the energy density of the turbulent field immediately
upstream of the shock. The most important finding is then that for Bohm
scaling of diffusion in the precursor with factor $\eta$, the change in
spectral index is invariably

$\Delta
s\lesssim\frac{2\,(s-1)^{2}\,\eta}{3\,\epsilon\,M_{\mathrm{A}}}={\color[rgb]{0,0,0}\frac{2\,(s-1)^{2}\,\eta}{3\,\epsilon}\sqrt{\frac{U_{\delta
B}}{U_{\mathrm{bulk}}}}}.$ (26)

Here $s$ is the cosmic-ray spectral index without steepening, and $\epsilon$
is the magnetic fraction of the turbulence energy density immediately upstream
of the shock. In Figure 1 we show this constraint as red exclusion area in a
display of spectral softening vs. $M_{\mathrm{A}}$.

The magnetic-field strength immediately upstream of the shock, and hence
$M_{\mathrm{A}}$, cannot be easily determined through observations.
Essentially all radiation is produced downstream of the shock, and there are
numerous processes operating at the shock that generate magnetic fields and
can cause a jump in $\mathrm{rms}(B)$ that is much larger than the compression
imposed at the shock (e.g. Giacalone & Jokipii, 2007; Fraschetti, 2013). The
true Alfvénic Mach number would then be larger than naively estimated.

The scale factor of the diffusion coefficient, $\eta$, could be large, but not
too large at very high energies, otherwise a young SNR could not accelerate
particles to that energy. Simulations suggest that $\eta\approx 4$, if only
the nonresonant mode is driven (Reville & Bell, 2013; Kobzar et al., 2017).
Waves driven by cosmic rays in one energy band also scatter cosmic rays with
different energies, albeit with somewhat modified efficiency (Reville et al.,
2008; Beresnyak & Li, 2014), which suggests that $\eta\approx 4$ is also a
reasonable estimate when one considers the entire spectrum of cosmic rays
accelerated at the shock. Modes driven by cosmic rays of the highest energies
may provide large-angle scattering for lower-energy cosmic rays, for which the
size of individual magnetic filaments is commensurate with the Larmor radius
in those filaments. Any other MHD mode, including those that are nonlinearly
driven by Bell’s mode, could lower $\eta$, because it may provide cosmic-ray
scattering and can reduce the Larmor radius, $r_{\mathrm{L}}$, through further
magnetic-field amplification. Likewise, MHD processes that generate large-
scale magnetic field in the precursor (e.g. Beresnyak et al., 2009; Bykov et
al., 2011; del Valle et al., 2016; Xu & Lazarian, 2017) would reduce $\kappa$
and hence $\Delta s$.

A remaining possibility is that the magnetic field in the precursor is not
amplified, in which case one might expect $\eta\approx B^{2}/(\delta B)^{2}$.
Inserting this into eq. 24 and using eq. 25 again yields $\Delta s$ of the
order of $M_{\mathrm{A}}^{-1}$. We conclude that $\eta$ can likely not
significantly increase the level of spectral softening.

Bell et al. (2019) based their estimate on the energy flux of turbulence and
cosmic rays through the shock, and the finding $\Delta s\propto
v_{\mathrm{sh}}$ is a consequence of the assumed saturation level $U_{\delta
B}\propto v_{\mathrm{sh}}$. Hadronic gamma-ray emission from SNRs shows soft
spectra mainly for older SNRs whose shocks have decelerated to about $0.01c$
or less, for which eq. 2 yields negligible steepening. Radio data are much
more abundant, but they reflect electron spectra. Technically,
ultrarelativistic electrons can also drive Bell’s mode, albeit with opposite
circular polarization (Bret, 2009). As cosmic-ray ions typically carry more
energy density than do the electrons, they would provide most of the amplified
field, and it is questionable whether any estimate of spectral steepening can
be directly applied to electrons and hence to radio synchrotron spectra.
Besides, the youngest known SNR in the Galaxy, G1.9+0.3, has a very fast
forward shock (Reynolds et al., 2008), but the radio emission most likely
comes from the reverse shock whose speed is considerably lower (Brose et al.,
2019).

Figure 1: Exclusion limits for the spectral softening, $\Delta s$, as a
function of the Alfvénic Mach number, $M_{\mathrm{A}}$, assuming $\eta=4$ and
$\epsilon=0.5$. The red area indicates a violation of eq. 26. The cyan and
blue areas are excluded by eq. 29 for
$v_{\mathrm{sh}}U_{\mathrm{cr}}/cU_{\mathrm{bulk}}$ set to $10^{-3}$ and
$4\cdot 10^{-3}$, respectively. The black solid line indicates the softening
that is imposed by the magnetic pressure.

It is instructive to assess the energy losses of cosmic-ray electrons. The
ratio of the energy-loss time of protons (Eq. 8) and the synchrotron loss time
of electrons at the same energy is independent of the magnetic-field strength.
The nonresonant streaming instability operates only upstream of the shock,
whereas the energy loss electrons incur in the compressed magnetic field
during the downstream half-cycle of shock acceleration is about ten times that
of the upstream half-cycle. Then, for parameters typical for young SNRs, a
shock speed $v_{\mathrm{sh}}=7000\,\mathrm{km\,s^{-1}}$,
$\omega_{\mathrm{p,p}}=1\,\mathrm{kHz}$, and
$N_{\mathrm{p}}=0.5\,\mathrm{cm^{-3}}$, the energy-loss time ratio exceeds
unity for 15-TeV particles and increases $\propto\Gamma_{\mathrm{cr}}^{2}$. At
50 TeV the synchrotron loss timescale of electrons would be 10% of the proton
loss time for driving Bell’s mode. Electron acceleration to 50 TeV or more
would imply that the acceleration timescale is still significantly shorter
than the synchrotron loss time, and so the energy loss of cosmic-ray protons
is negligible at this energy. In other words, if intense nonthermal multi-keV
emission is produced, the energy loss, and hence the spectral steepening, is
very small for hadronic cosmic rays that produce TeV-band gamma-ray emission.

Instead of a fixed saturation level, we consider the energy transfer
throughout the cosmic-ray precursor, which accounts for the limited available
time and includes energy that nonlinearly passes to heating or other types of
turbulence, but requires the assumption that the linear growth rate of the
non-resonant mode is also valid in the nonlinear stage. A consistency check is
provided by mandating that in the steady state the wave-driving power per unit
shock area (eq. 21) must be at least as large as the escape flux through the
shock, $v_{\mathrm{sh}}U_{\delta B}/\epsilon$. This yields an implicit lower
limit on the cosmic-ray density that with eq. 5 corresponds to requiring at
least one half exponential growth cycle for the nonresonant mode,

$U_{\mathrm{cr}}\geq\frac{6}{\eta}\,\frac{U_{\mathrm{bulk}}}{M_{\mathrm{A}}}\quad\Rightarrow\
N_{\mathrm{exp}}\gtrsim 0.5.$ (27)

Evidently that condition is met for any significant growth of the mode. As the
mode has to grow from small fluctuations to a level $\delta B\gg B_{0}$ , one
would instead need around ten exponential growth cycles. Interestingly, we
recover the result of Bell et al. (2019), $\Delta s|_{\mathrm{Bell}}\approx
4\,U_{\delta B}/(\epsilon U_{\mathrm{cr}})$, if we assume a cosmic-ray density
that equals the right-hand side eq. 27. In that case $N_{\mathrm{exp}}\simeq
0.5$, meaning there is no time for the mode to grow, and so the steady-state
level $U_{\delta B}/U_{\mathrm{cr}}\approx v_{\mathrm{sh}}/(2c)$ cannot be
reached.

Our calculation of the number of exponential growth times in eq. 5 implicitly
assumes that the mean diffusion coefficient, $\kappa$, in the cosmic-ray
precursor is not much different from that immediately upstream of the shock,
at least for cosmic rays at energies well below the cut-off. We noticed in
section 2.2 that this must by roughly true, otherwise most upstream cosmic
rays would reside very far from the shock and effectively escape, which is a
far more serious loss process than wave driving.

The spectral modification implied by Eq. 26 does not explicitely depend on the
energy density of cosmic rays, $U_{\mathrm{cr}}$. The spectral softening
reflects the per-particle rates of energy gain and escape from the shock. The
energy transfer rate from cosmic rays to turbulence scales linearly with the
cosmic-ray density, but normalized to the individual cosmic ray it is
independent. Hence the independence of $\Delta s$ on the cosmic-ray density.
We calculated the level of spectral softening ignoring the condition
$kr_{\mathrm{L}}\gg 1$, under which the non-resonant mode can be driven. For
$\delta B\gg B_{0}$ this relation can be rephrased as

$\frac{v_{\mathrm{sh}}}{4c}\frac{U_{\mathrm{cr}}}{U_{\delta B}}\gg 1.$ (28)

This condition requires that the left-hand side be much larger than unity, but
how much larger? Let us conservatively suppose that is it larger than or equal
to two. Inserting that into eq. 26 we find

${\color[rgb]{0,0,0}\Delta
s\lesssim\frac{(s-1)\eta}{3\sqrt{2}\,\epsilon}\sqrt{\frac{v_{\mathrm{sh}}}{c}\frac{U_{\mathrm{cr}}}{U_{\mathrm{bulk}}}}}.$
(29)

This condition is also displayed in Figure 1 for two rather high values of
$v_{\mathrm{sh}}$ and $U_{\mathrm{cr}}$. It appears as a constant limit at all
$M_{A}$, although some values of $M_{A}$ may not be reached, at least not by
magnetic-field amplification through Bell’s mode.

Given that the standard jump condition increases $\sqrt{U_{\delta
B}/U_{\mathrm{bulk}}}$ in the downstream region to about 8 times its value in
the upstream region, a strong magnetic field, or small $M_{\mathrm{A}}$, will
by itself modify the shock (e.g. Vainio & Schlickeiser, 1999; Lerche et al.,
2000; Caprioli et al., 2009) and hence soften the particle spectrum in a
similar way as does the energy loss by driving turbulence. The solid black
line in Figure 1 indicates the magnitude of this effect for a quasi-
perpendicular shock with dynamically relevant field amplitude of
$\sqrt{2/3}\,\delta B$. There is a marginal dependence of the curve on the
shock-parallel magnetic-field component that we ignore here. Note also that a
turbulently oriented magnetic field will for $M_{\mathrm{A}}\lesssim 20$
induce vorticity at the shock that can drive a turbulent dynamo in the
downstream region and further amplify the field there.

In conclusion, we find and show in Figure 1 that even for very efficient
cosmic-ray acceleration, for which $\eta\approx 4$, and the highest magnetic-
field amplification that is allowed for Bell’s mode, the spectral softening
appears to be moderate, $\Delta s\lesssim 0.1$, and it is negligible for
standard SNR parameters. As we explicitely allow for spatial variation in the
cosmic-ray precursor, the shock speed in Eq. 29, $v_{\mathrm{sh}}$, is that of
the thermal subshock, not that measured in the far-upstream frame. This
statement is based on the energy transfer that can be accomplished within the
time available, and it does not assume that a certain saturation level of the
wave energy density is reached.

## References

  * Amato & Blasi (2006) Amato, E., & Blasi, P. 2006, MNRAS, 371, 1251, doi: 10.1111/j.1365-2966.2006.10739.x
  * Bell (1978) Bell, A. R. 1978, MNRAS, 182, 147, doi: 10.1093/mnras/182.2.147
  * Bell (2004) —. 2004, MNRAS, 353, 550, doi: 10.1111/j.1365-2966.2004.08097.x
  * Bell et al. (2019) Bell, A. R., Matthews, J. H., & Blundell, K. M. 2019, MNRAS, 488, 2466, doi: 10.1093/mnras/stz1805
  * Bell et al. (2013) Bell, A. R., Schure, K. M., Reville, B., & Giacinti, G. 2013, MNRAS, 431, 415, doi: 10.1093/mnras/stt179
  * Beresnyak et al. (2009) Beresnyak, A., Jones, T. W., & Lazarian, A. 2009, ApJ, 707, 1541, doi: 10.1088/0004-637X/707/2/1541
  * Beresnyak & Li (2014) Beresnyak, A., & Li, H. 2014, ApJ, 788, 107, doi: 10.1088/0004-637X/788/2/107
  * Bret (2009) Bret, A. 2009, ApJ, 699, 990, doi: 10.1088/0004-637X/699/2/990
  * Brose et al. (2020) Brose, R., Pohl, M., Sushch, I., Petruk, O., & Kuzyo, T. 2020, A&A, 634, A59, doi: 10.1051/0004-6361/201936567
  * Brose et al. (2019) Brose, R., Sushch, I., Pohl, M., et al. 2019, A&A, 627, A166, doi: 10.1051/0004-6361/201834430
  * Bykov et al. (2011) Bykov, A. M., Osipov, S. M., & Ellison, D. C. 2011, MNRAS, 410, 39, doi: 10.1111/j.1365-2966.2010.17421.x
  * Caprioli et al. (2009) Caprioli, D., Blasi, P., Amato, E., & Vietri, M. 2009, MNRAS, 395, 895, doi: 10.1111/j.1365-2966.2009.14570.x
  * del Valle et al. (2016) del Valle, M. V., Lazarian, A., & Santos-Lima, R. 2016, MNRAS, 458, 1645, doi: 10.1093/mnras/stw340
  * Drury (1991) Drury, L. O. 1991, MNRAS, 251, 340, doi: 10.1093/mnras/251.2.340
  * Fraschetti (2013) Fraschetti, F. 2013, ApJ, 770, 84, doi: 10.1088/0004-637X/770/2/84
  * Gargaté et al. (2010) Gargaté, L., Fonseca, R. A., Niemiec, J., et al. 2010, ApJ, 711, L127, doi: 10.1088/2041-8205/711/2/L127
  * Giacalone & Jokipii (2007) Giacalone, J., & Jokipii, J. R. 2007, ApJ, 663, L41, doi: 10.1086/519994
  * Jansson & Farrar (2012) Jansson, R., & Farrar, G. R. 2012, ApJ, 761, L11, doi: 10.1088/2041-8205/761/1/L11
  * Kobzar et al. (2017) Kobzar, O., Niemiec, J., Pohl, M., & Bohdan, A. 2017, MNRAS, 469, 4985, doi: 10.1093/mnras/stx1201
  * Lerche et al. (2000) Lerche, I., Pohl, M., & Schlickeiser, R. 2000, Journal of Plasma Physics, 64, 459, doi: 10.1017/S0022377800008801
  * Luo & Melrose (2009) Luo, Q., & Melrose, D. 2009, MNRAS, 397, 1402, doi: 10.1111/j.1365-2966.2009.14872.x
  * Niemiec et al. (2008) Niemiec, J., Pohl, M., Stroman, T., & Nishikawa, K.-I. 2008, ApJ, 684, 1174, doi: 10.1086/590054
  * Ohira et al. (2010) Ohira, Y., Murase, K., & Yamazaki, R. 2010, A&A, 513, A17, doi: 10.1051/0004-6361/200913495
  * Reville & Bell (2013) Reville, B., & Bell, A. R. 2013, MNRAS, 430, 2873, doi: 10.1093/mnras/stt100
  * Reville et al. (2009) Reville, B., Kirk, J. G., & Duffy, P. 2009, ApJ, 694, 951, doi: 10.1088/0004-637X/694/2/951
  * Reville et al. (2007) Reville, B., Kirk, J. G., Duffy, P., & O’Sullivan, S. 2007, A&A, 475, 435, doi: 10.1051/0004-6361:20078336
  * Reville et al. (2008) Reville, B., O’Sullivan, S., Duffy, P., & Kirk, J. G. 2008, MNRAS, 386, 509, doi: 10.1111/j.1365-2966.2008.13059.x
  * Reynolds et al. (2008) Reynolds, S. P., Borkowski, K. J., Green, D. A., et al. 2008, ApJ, 680, L41, doi: 10.1086/589570
  * Riquelme & Spitkovsky (2009) Riquelme, M. A., & Spitkovsky, A. 2009, ApJ, 694, 626, doi: 10.1088/0004-637X/694/1/626
  * Romanova & Owocki (2015) Romanova, M. M., & Owocki, S. P. 2015, Space Sci. Rev., 191, 339, doi: 10.1007/s11214-015-0200-9
  * Stroman et al. (2009) Stroman, T., Pohl, M., & Niemiec, J. 2009, ApJ, 706, 38, doi: 10.1088/0004-637X/706/1/38
  * Vainio & Schlickeiser (1999) Vainio, R., & Schlickeiser, R. 1999, A&A, 343, 303
  * Voelk & Forman (1982) Voelk, H. J., & Forman, M. 1982, ApJ, 253, 188, doi: 10.1086/159623
  * Weidl et al. (2019) Weidl, M. S., Winske, D., & Niemann, C. 2019, ApJ, 873, 57, doi: 10.3847/1538-4357/ab0462
  * Xu & Lazarian (2017) Xu, S., & Lazarian, A. 2017, ApJ, 850, 126, doi: 10.3847/1538-4357/aa956b


Manganese in the West Atlantic Ocean in context of the first global
ocean circulation model of manganese
[1,6]Marcovan Hulten
[4,5]Heinde Baar

[1]Laboratoire des Sciences du Climat et de l'Environnement (LSCE), IPSL,
CEA–Orme des Merisiers, 91191 Gif-sur-Yvette, France
[2]Department of Chemistry, NIWA/University of Otago Research Centre for
Oceanography, Dunedin 9054, New Zealand
[3]Department of Ocean Sciences & Institute of Marine Sciences,
University of California Santa Cruz, CA 95064, USA
[4]NIOZ Royal Netherlands Institute for Sea Research, Department of Ocean
Systems, and Utrecht University, P.O. Box 59, 1790 AB Den Burg, Texel, the
[5]University of Groningen (RUG),
Postbus 72, 9700 AB Groningen, the Netherlands
[6]Royal Netherlands Meteorological Institute (KNMI),
Utrechtseweg 297, 3731 GA De Bilt, the Netherlands
[7]University of Liverpool,
4 Brownlow Street, Liverpool L69 3GP, UK

Manganese in the Atlantic Ocean
M.M.P. van Hulten et al.
M. M. P. van Hulten ()

Dissolved manganese (Mn) is a biologically essential element.
Moreover, its oxidised form is involved in removing itself and several other trace
elements from ocean waters.
Here we report the longest thus far 17 500 km length full-depth ocean section
of dissolved Mn in the West Atlantic Ocean, comprising 1320 data values of high
This is the GA02 transect that is part of the Geotraces
programme, which aims to understand trace element distributions.
The goal of this study is to combine these new observations with a new,
state-of-the-art, modelling to give a first assessment of the main sources and
redistribution of Mn throughout the ocean.
To this end, we simulate the distribution of dissolved Mn using a
global-scale circulation model.
This first model includes simple parameterisations to account
for the sources, processes and sinks of Mn in the ocean.
Oxidation and (photo)reduction, aggregation, settling, as well as biological
uptake and remineralisation by plankton, are included in the model.
Our model provides, together with the observations, the following insights:

* The high surface concentrations of manganese are caused by the
combination of photoreduction and sources to the upper ocean.
The most important sources are sediments, dust, and, more locally,

* Observations and model simulations suggest that surface Mn in the Atlantic
Ocean moves downwards into the southward flowing North Atlantic Deep Water
(NADW), but because of strong removal rates there is no elevated concentration
of Mn visible any more in the NADW south of 40 N.

* The model predicts lower dissolved Mn in surface waters of the
Pacific Ocean than the observed concentrations.
The intense Oxygen Minimum Zone (OMZ) in subsurface waters is deemed to be a
major source of dissolved Mn also mixing upwards into surface waters, but the
OMZ is not well represented by the model.
Improved high resolution simulation of the OMZ may solve this problem.

* There is a mainly homogeneous background concentration of dissolved Mn of
about 0.100.15 throughout most of the deep ocean.
The model reproduces this by means of a threshold on particulate manganese
oxides of 25, suggesting that a minimal concentration of
particulate Mn is needed before aggregation and removal become efficient.

* The observed distinct hydrothermal signals are produced by assuming both a
strong source and a strong removal of Mn near hydrothermal vents.

Dissolved manganese (Mn_diss) is taken up by phytoplankton,
because Mn is crucial for photosynthesis and other biological functions
Furthermore, its oxidised form
(Mn_ox) plays an important role in the removal of several
other trace metals from seawater <cit.>.
While in the open ocean manganese (Mn) exists in small concentrations, it is the
twelfth most plentiful element in the Earth's crust <cit.>.
In seawater, Mn occurs in many forms, among which the bioavailable dissolved
After phytoplankton death, incorporated Mn sinks downwards
together with the dead material, but most of the organic material is
remineralised before reaching the sea floor <cit.>,
releasing Mn back to the water.

Another important mechanism of storing Mn in particles, besides biological
incorporation, is the removal of dissolved Mn via larger colloids on which the Mn
is oxidised to insoluble Mn(IV) (and possibly the other Mn(III) oxidation
state), and the subsequent aggregation by particulate matter.
occurs everywhere in the ocean where oxygen is available.
This process can be
strongly accelerated by Mn(II)-oxidising microorganisms, primarily bacteria and
fungi <cit.>, but we do not understand the role of
these organisms quite well <cit.>.
The reverse process is the reduction of Mn oxides to bioavailable
dissolved Mn(II), i.e. Mn$^{2+}\text{({\it aq})}$.
The full oxidation/reduction (redox) equilibrium reaction, in its most simple
form, is given by <cit.>:

Mn^2+ + 1/2 O_2 + H_2O
MnO_2 + 2 H^+   .

Reduction is significantly faster under the influence of sunlight.
It is hence referred to as photoreduction when irradiance is the major
The relative rate of reduction compared to oxidation
is important for Mn(II) availability.
An overall net higher oxidation rate implies more particle formation, hence more Mn export.

Manganese enters the open ocean through lithogenic dust deposition <cit.> and
lateral advection from reducing sediments <cit.>.
Sediments along the relatively shallow (<1000 m depth) ocean margins tend to
receive more organic deposition, hence by enhanced microbial decomposition are
more anoxic and a stronger source of reduced Mn than deep sea sediments.
Surface [Mn_diss] is especially high in the central Atlantic Ocean and up
to at least 30 N
because of high dust input from the Sahara in combination with photoreduction
[Mn_diss] is high in the northern Indian Ocean <cit.>.

Dissolved Mn diffuses out of oxygen deprived sediments, because sediment microorganisms reduce
Mn_ox if there is no more oxygen or nitrate left
As long as there is oxygen in the sediment, the organic carbon is remineralised by
using this oxygen as an electron acceptor <cit.>.

Rivers are another source of Mn to the ocean <cit.>.
While much of the fluvial Mn is removed within the estuaries into the sediments by scavenging and
aggregation, a large part may finally be transported to the ocean by diffusion
from and resuspension of the sediments <cit.>.
Typically, the smaller sediment particles (0.54)
have a high Mn content, and, because of their small size, are able to reach
the open ocean <cit.>.
Manganese may also flux into the ocean by melting sea ice <cit.>.
Finally, overwhelming evidence is found of manganese fluxing out of hydrothermal vents

Downward fluxes of settling particles that have been collected in sediment traps show a
strong correlation between lithogenic particles and authigenic Mn
Therefore, lithogenic particles are likely to play a
significant role in the removal (oxidation, scavenging and aggregation) of Mn
Here the ballast effect of lithogenic particles, which typically have a density
of about twice of that of
seawater, is likely playing a major role for rapid settling of
agglomerates of biogenic and lithogenic particles including Mn-oxides coatings
This is consistent with the fact that most suspended particles are small (less
than 2), but the larger aggregates are deemed to be the
significant contributors to the vertical flux <cit.>.
The complete process may be more complicated than described above
<cit.>, e.g. because Mn binds to dissolved ligands such that more of it may
stay in solution <cit.>.

Manganese oxides are important scavengers of other trace metals like iron,
cobalt, nickel and zinc <cit.>, as well as
insoluble radionuclides such as thorium and protactinium

Therefore, Mn availability does not only directly impact primary production but
may also play a role in removing other elements from the surface ocean.
These elements include biologically essential trace metals as well as
many more trace elements.

Published observational studies, like the results of <cit.>,
show low, relatively constant Mn_diss concentrations of around
0.15 away from the boundaries.

Towards the surface there is a sharp increase of [Mn_diss].
Other notable elevations are near oceanic ridges.

These elevations of [Mn_diss] at mid-depths can be ascribed to
hydrothermal activity near those regions.
The hydrothermal plumes typically extend in the order of 1000
<cit.>, or reach even up to
3000 <cit.>.

What makes the distribution of Mn_diss relatively homogeneous in
the interior of the ocean?
How can this be reconciled with localised features where [Mn_diss] is
very elevated compared to that stable “background concentration”?

Until now only local model simulations of the Mn ocean cycle have been
performed that focus on the processes most relevant for the respective
In other words, the different features of the Mn_diss distribution
have not been brought together in a unified model.

The case of hydrothermal activity has been studied by
<cit.> who modelled Mn in the deep ocean near hydrothermal
They included four Mn tracers, namely, a dissolved form, small particles
associated with bacteria, larger aggregate particles, and one in
sediments as the model includes benthic fluxes.
They found that “more than 80 % of the hydrothermal Mn is deposited
within several hundred kilometres of the ridge crest though dissolved Mn
concentrations beyond that distance exceed background levels by many
This illustrates a high Mn input and removal from hydrothermal vents, as well
as the already mentioned large plume extend.

The Mn_diss distribution in the North Pacific Ocean has
been modelled by <cit.>, in a 1-D vertical model neglecting
horizontal transport.
Their oxidation model depends on the oxygen concentration and the hydroxide activity.
Their region of interest was on the upper and intermediate depth ocean,
and their goal was to reproduce the [Mn_diss] maximum in
the Oxygen Minimum Zone (OMZ).
Besides the North Pacific, the OMZ is also present in other basins, amongst which
the northwestern Indian Ocean <cit.>.
<cit.> found that the combination of remineralisation rates and
decreased oxidation in the OMZ explained the [Mn_diss]
subsurface maximum in their profiles, whereas a flux from the continental margin
sediments appeared not responsible.
In the euphotic zone Mn_diss was incorporated in
phytoplankton, while remineralised Mn_diss in the aphotic
(and disphotic) zone was lost by oxidation and scavenging.

While these modelling efforts are useful for their purposes, no studies
exist in which the global ocean Mn_diss distribution is modelled.
To arrive at a more integrated understanding of the behaviour of pelagic Mn,
we include a Mn model in a global ocean general circulation model.
Specifically, we will test a simple mechanism that should be able to give
insight into the apparently stable Mn_diss concentration, and its
contrast to strong Mn sources.
This is the first time that a global ocean model for manganese has been written and
It is a basic model that should give a starting point for further studies.
At the moment there is no mechanistic evidence for typical
uptake-remineralisation processes as is the case for, e.g., iron.
However, as there is of course uptake of Mn by phytoplankton, and sometimes it
can even be a limiting nutrient <cit.>, we include a biological
cycle of Mn in the model.
While we perform a model simulation on a global scale, we give more
attention to the Mn_diss distribution in the Atlantic Ocean.
We will compare our simulations in most detail, and quantitatively, with a
here reported highly accurate dataset from the Geotraces programme,
namely the GA02 section in the West Atlantic Ocean.
Here accuracy refers to the proximity of the measurements to the true
values, as confirmed by (i) the compliance of reference samples with
international consensus values (Table 1), and (ii) the agreement at the crossover
station (Fig. <ref>) with the US section GA03, and (iii) the
agreement of the two independent methods FIA and ICPMS for the our data of
section GA02 (Fig. <ref>).
In addition, the North Atlantic GA03 <cit.> and the Zero Meridian
Southern Ocean GIPY5 <cit.> transects will be used for
further detailed visual comparison.
Furthermore, in those regions we can study important properties of the ocean
geochemistry of manganese and the interaction with circulation, among which the
Atlantic overturning circulation and hydrothermal activity.

In this study our goal is to assess the fundamental processes that are the most
important to accurately simulate [Mn_diss], including
the aforementioned properties of the dissolved Mn distribution.
In this case, accuracy refers to the proximity of the model to the
For clarity we will always qualify if this relates to the closeness of the model
simulation to the observations instead of the proximity of the measurements to
the true values.
To this end, we will first introduce our Mn model with its processes, sources and sinks.
We will show the results of a reference simulation, which will be
compared with recent high-accuracy observations.
Also several sensitivity simulations will be presented, studying the effects of
the biological cycle, the intense nature of hydrothermal vents and the strong
removal of manganese from the seawater.

§ METHODS

§.§ Observations in the West Atlantic Ocean

§.§.§ Sample collection

For the determination of trace metal concentrations,
samples were collected along the Geotraces Atlantic
Meridional GA02 transect of the Netherlands (Fig. <ref>).
Sampling was done with an all-titanium ultraclean CTD sampling
system for trace metals <cit.> with novel PVDF samplers
Immediately upon recovery, the complete titanium frame with its 24 PVDF samplers
was placed inside a clean room environment where the sub-samples
for trace metal analysis were collected. The water was filtered from the PVDF
samplers over a 0.2 filter cartridge (Sartobran-300,
Sartorius) under pressure (1.5 atm) of (inline prefiltered) nitrogen gas.
Sub-samples for dissolved metals were taken in cleaned <cit.> LDPE sample bottles. All sample bottles were rinsed five
times with the sample seawater. Seawater samples were acidified with HCl to a
concentration of 0.024 M HCl which results in a pH of 1.7 to 1.8 with
Baseline® Hydrochloric Acid (Seastar Chemicals Inc.).

§.§.§ Analysis of dissolved Mn

Analyses of dissolved manganese were performed shipboard with the method
developed by <cit.> with some slight modifications in the preparation
and brands of the chemicals used.
Notably, samples were buffered in-line
with an ammonium borate sample buffer to a pH of $8.5 \pm 0.2$.
This buffer was produced by dissolving
30.9 g of boric acid (Suprapure, Merck) in 1 L MQ water (Millipore Milli-Q)
deionised water $R > \SI{18.2}{\mega\ohm\per\centi\metre}$ and adjusting the pH
to 9.4 with ammonium hydroxide (Suprapure, Merck).

The buffered sample was pre-concentrated during 150 s on a Toyopearl AF-Chelate
650M (TosoHaas, Germany) column. Hereafter the column was rinsed for 60 s with
MQ water to remove interfering salts. The Mn was
subsequently eluted from the column for 200 seconds with a solution of 0.1 M
three times quartz distilled formic acid (reagent grade, Merck) containing
0.1 M hydrogen peroxide (Suprapure, Merck) and 12 mM ammonium hydroxide
(Suprapure, Merck). The pH of this carrier solution was adjusted to $2.9 \pm
0.05$. The eluate with the dissolved Mn passed a second column of immobilised
8-hydroxyquinoline <cit.> to remove interfering iron ions in the
carrier solution <cit.>. Hereafter the carrier mixed
with 0.7 M ammonium hydroxide (Suprapure Merck) and a luminol solution. The
latter luminol solution was made by diluting 600 luminol
stock solution and 10 TETA (triethylenetetramine, Merck) in
1 MQ. The luminol stock solution was made by diluting 270 mg luminol
(3-aminophtalhydrazide, Aldrich) and 500 mg potassium carbonate in 15 ml MQ.
The resulting mixture of carrier solution, ammonium hydroxide and luminol
solution had a pH of $10.2 \pm 0.05$ and entered a 3 m length mixing coil
placed in a water bath of 25 C. Hereafter the chemiluminescence was
detected with a Hamamatsu HC135 Photon counter.  Concentrations of dissolved Mn
were calculated in nanomole per litre (nM) from the photon emission peak height
of triplicate measurements.

The system was calibrated using standard additions from a 5000 nM Mn stock
solution (Fluka) to filtered acidified seawater of low Mn concentration that was
collected in the sampling region. A five-point calibration line (0, 0.1, 0.2, 0.6
and 1.2 nM standard additions) and blank determination were made every day. The
three lowest points (0, 0.1 and 0.2 nM) of the calibration line were measured
in triplicate and the two highest points (0.6 and 1.2 nM) in duplicate in order
to add more weight to the lower part of the calibration line.  The blank was
determined by measuring acidified MQ which was below the detection limit and
subsequently no blank was substracted. The limit of detection defined as three
times the standard deviation of the lowest value observed was $< 0.01$ nM. The
flow injection system was rinsed every day with a 0.5 M HCl solution.

An internal reference sample was measured in triplicate every day. This was a
sub-sample of a 25 volume of filtered seawater that was taken at the
beginning of Leg 1 (also used during Leg 2, i.e. Iceland to the equator) and Leg 3
(Punta Arenas, Chile, to the equator). The relative standard deviation (i.e.the precision) of this replicate analysis seawater sample that was analysed 40
times on different days in triplicate was 2.57 % (Leg 1 and 2) and 1.21 %
for 17 analyses during Leg 3. The relative standard deviation on single days was
on average 1.37 % and the absolute values were 0.45 and 0.61 nM for the first
two legs and Leg 3, respectively.

Mn_diss ()     Consensus     <cit.>

SAFe D2                                 $0.35 \pm 0.05  $       $0.33\pm 0.01 \;(n=24)$

Geotraces S                    $1.46 \pm 0.14  $       $1.47\pm 0.03 \;(n=10)$

Geotraces D                    $0.21 \pm 0.03  $       $0.18\pm 0.01 \;(n= 5)$

Compliance with the international Geotraces reference samples program.
Left column are the international consensus values. Right column are the values
reported here as part of the GA02 dataset.

As external comparison, the international reference samples collected on the
Geotraces Intercalibration Cruise (<www.geotraces.org>) as well as
from the SAFe cruise <cit.> were analysed for Mn.

At different locations, large volumes of seawater were sampled of which
subsamples were analysed by various laboratories, resulting in different
reference values.
These labs include the Royal Netherlands Institute for Sea Research (NIOZ), who
also analysed the GA02 transect data.
An independent referee removed outliers and averaged the reference values,
resulting for each of the sample locations in a consensus value,
which is here considered as the true value.
These consensus values are listed in Table <ref>, together with
the values determined by NIOZ.
The distributions of the measurements from NIOZ lie within the precision of the
consensus values.

[Mn_diss] at the 16 June 2010 (NL) and the
19 November 2011 (US) occupations of the crossover site of the Bermuda Atlantic
Time Series (BATS) station.
There is good agreement between the two sets of samples, as well as the two
different analytical methods, the shipboard FIA by NIOZ for the NL2010 samples,
and the lab ICP-MS at UCSC for both the NL2010 and the US2011 samples.
Each lab had its own independently prepared lab standards, further confirming
the overall accuracy.
Analyses by Rob Middag at NIOZ in 2010 and in 2011 at UCSC.

The vertical profiles of the BATS station for the Netherlands sampling at
13 June 2010 and the US sampling at 19 November 2011 are presented in
Fig. <ref>.
For both occupations the samples were analysed by ICP-MS by Rob Middag at UCSC.
Moreover, the same analyst had done shipboard FIA analyses during the 2010 cruise.
There was no statistical difference within analytical uncertainty between the
three datasets in the deep waters <cit.>.
Small differences within the upper $\sim$500 m between the Netherlands sampling
and the US sampling are deemed to be real due to the seasonal variability over
the more than one year time difference of sampling.
Furthermore, these Mn profiles are also consistent with the profiles determined
by <cit.>.

The correlation between the two methods of analysis for the
determination of [Mn_diss] at GA02 (shipboard and laboratory
mass spectrometry): (a) at Bermuda; (b) all 55 West Atlantic stations.

Besides the shipboard Flow Injection Analysis (FIA), we also analysed samples
using mass spectrometer analysis.
Figure <ref> shows the correlation between the two methods of
analysis for the determination of [Mn_diss], shipboard and
laboratory measurements.
There is a very good agreement between the shipboard and mass spectrometer
analyses, which strongly suggests a high observational accuracy.
For the GIPY5 transect in the Southern Ocean, we analysed
[Mn_diss] only through flow injection (FIA).
Since we want to plot GA02 together with the GIPY5 in this study, we use the
observational data obtained from the shipboard FIA for the comparison with the
model simulations.

§.§ Model description

In order to simulate the three-dimensional (3-D) distribution of dissolved Mn,
we use the general circulation model Océan PArallélisé (OPA)
that is part of NEMO, a framework for ocean models <cit.>.
We use the ORCA2-LIM configuration of NEMO.
The spatial resolution is 2$^{\circ}$ by 2$^{\circ}\,\cos(\phi)$ (where $\phi$
is the latitude) with an increased meridional resolution to 0.5$^{\circ}$ in the
equatorial domain <cit.>.
The model has 30 vertical layers, with an increased vertical thickness from
10 at the surface to 500 at 5000 depth.
Representation of the topography is based on the partial-step thickness
Lateral mixing along isopycnal surfaces is performed both on tracers and
momentum <cit.>.
The parameterisation of <cit.> is applied from 10 poleward to
represent the effects of non-resolved mesoscale eddies.
Vertical mixing is modelled using the turbulent kinetic energy (TKE) scheme of
<cit.>, as modified by <cit.>.
The fluid dynamics used to drive our model is identical to that used in

The model contains two tracers of Mn, referred to as dissolved
(Mn_diss) and oxidised (Mn_ox) manganese.
These tracers are driven by the equations set out in this section, as well as
the velocity fields obtained from OPA.
Instead of calculating the dynamical variables of the model (velocity and
mixing), we run it off-line, using a climatology with a resolution of five days.

Model scheme:
biology in green; redox and scavenging in blue; circulation and mixing
(dynamics) in grey; and the sources dust, rivers, hydrothermal and sediment
are in light brown, light blue, dark brown and red, respectively.
Mn_diss is the dissolved and Mn_ox the oxidised Mn.

Figure <ref> presents the conceptual scheme of our
manganese model.
The internal processes include biological uptake and remineralisation (green
arrows), reduction, oxidation, aggregation and burial (blue arrows in the
figure), and are described in the following subsections.
Manganese sources from rivers, the atmosphere, sediments and hydrothermal vents
are presented as arrows at the top and bottom of the figure.
These four Mn sources are presented in Fig. <ref>, and
Table <ref> lists the absolute
contributions to the different basins by each of these sources, as well as
the relative contribution of every source to the world ocean.
The model parameters are summarised in Table <ref>.
In the following subsections we will describe how the different sources and
processes are included in the model.

Sources of Mn to the ocean: effective Mn_diss input flux $\Phi$.
Three-dimensional fields are vertically integrated, such that dimensions
are molar fluxes.
Ranges vary between the different sources.

Basin               Dust        Rivers      Sediment    Hydrothermal

Atlantic Ocean      2200          127            924          13917

Pacific Ocean       1183           94           1237          59846

Indian Ocean        1506           25            442          18653

Southern Ocean        14            0            206           7601

Arctic Ocean          23           19            463           2269

Mediterranean Sea     673          12             91              0

total amount (Mmol yr$^{-1}$)
5598          277           3363          102286

relative amount (%)     5.0      0.2           3.0           91.7

Absolute amount of effective annual input by an imposed Mn flux from
sediments, rivers and hydrothermal vents into each basin
The Southern Ocean is defined as the ocean south of 58.7 S.
The line “total amount” denotes how much world-wide Mn is added to the
ocean due to a specific flux; “relative amount” is normalised to the
total Mn input flux.

Parameter                       Symbol                        Value used                    Known range           Reference

Mass fraction of Mn in dust     $f_\text{Mn,dust}$            880 ppm                696880 ppm     <cit.>;

Dust Mn solubility              $\alpha$                       40 %                 1070 %        <cit.>

Sediment source Mn/Fe ratio     $r_\text{Mn:Fe,sed}$            0.2                   uncertain                 <cit.>;

River source Mn/Fe ratio        $r_\text{Mn:Fe,riv}$            0.214                 uncertain                 <cit.>

Hydrothermal Mn/^3He ratio     $r_{\chem{Mn:^3\!He},\text{hydro}}$     0.10e9     uncertain            -

Settling speed of Mn_ox     $w_\text{ox}$       1–10     0.9–1.4     <cit.>

Oxidation rate constant         $k_\text{ox}$                 0.341e-3      uncertain                               <cit.>;

Photoreduction rate constant    $k_\text{red,light}$          98e-3         50–150e-3        <cit.>

Aphotic reduction rate          $k_\text{red,dark}$           1.70e-3       0.98–14.3e-3     <cit.>;

Aggregation threshold           $\mathcal{X}_\text{thr}$      25          hypothetical                            -

Mn/P incorporation ratio        $R_\text{Mn:P}$               0.36e-3                 0.2e-31.5e-3              <cit.>;

Mn model parameters for the Reference simulation.
The settling velocity $w_\text{ox}$ is given by Eq. <ref>.
In LowHydro the hydrothermal flux and the maximum settling speed are
both reduced by a factor of 10 ($r_{\chem{Mn:^3\!He}\text{,hydro}} =
\num{0.01e9}$; and $w_\text{ox} \equiv \SI{1}{\metre\per\day}$ or,
equivalently, $\mathcal{X}_\text{thr} \rightarrow \infty$).
In the sensitivity simulation NoThreshold the aggregation threshold
$\mathcal{X}_\text{thr}$ is set to zero.

§.§.§ Atmospheric source

Manganese is added to the pool of Mn_diss in the upper model layer, according to

\begin{equation}
\frac{\partial \mathcal{D}}{\partial t}\Big|_\text{dust}
= \frac{\alpha \cdot f_\text{Mn,dust}}{m \cdot \Delta z_1}
\cdot \Phi_\text{dust} \,,
\end{equation}

where $\mathcal{D}$ is the dissolved Mn concentration, $\alpha$ is
the solubility of Mn in dust, $f_\text{Mn,dust}$ is the mass fraction of Mn
in dust, $m$ is the molar mass of Mn, and $\Delta z_1 =
\SI{10}{\metre}$ is the upper model layer thickness.
The lithogenic dust deposition flux, $\Phi_\text{dust}$, is derived from
the Interaction with Chemistry and Aerosols (INCA) model <cit.>.
Here we use a 12 month climatology of INCA's output as a forcing.

The average mass fraction of Mn in the Earth's upper
crust is 527 ppm <cit.>.
However, the fraction measured in Saharan dust is 880 ppm
<cit.>, consistent with <cit.> and
Since most of the dust deposited on the Atlantic Ocean originates from
the Sahara and our focus is the Atlantic Ocean, the value of 880 ppm is
used for $f_\text{Mn,dust}$.

The solubility of Mn from dust is uncertain and relatively high compared
to most other trace metals.
Here $\alpha = 40$ % of the Mn in dust is assumed to
dissolve, largely consistent with the values reported by
Several studies report even higher values ($> 50$ %), which are,
however, mainly from anthropogenic or otherwise processed dust, while
the lower reported values ($< 50$ %) are from natural dust, mainly of
Saharan origin.

Figure <ref>a presents the average Mn dissolution flux
of the 12 months climatology.
Globally this is 5.6 of Mn.

§.§.§ River source

The manganese river source is modelled analogously to iron, which is part of the
biogeochemical model Pisces <cit.>.
This means that our manganese influx is proportional to the total dissolved
(organic and inorganic) carbon flux, just like for iron <cit.>.
Hence, the modelled concentration change of Mn_diss caused by river input is
given by

\begin{equation}
\frac{\partial\mathcal{D}}{\partial t}\Big|_\text{rivers} =
r_\text{Mn:Fe,river} \cdot \frac{\partial[\chem{Fe}_\text{diss}]}%
{\partial t}\Big|_\text{rivers}\,,
\label{eqn:mang:rivers}
\end{equation}

where $[\chem{Fe}_\text{diss}]$ is the dissolved iron concentration, and
$r_\text{Mn:Fe,river}$ is the effective manganese/iron
flux ratio, here set to 0.214, based on river dissolved concentrations
For comparison, the Mn/Fe flux ratio is much higher than the crustal ratio,
which is around 0.02 <cit.>, but probably underestimated as we do
not consider external sources of particulate manganese.
The effective Mn_diss input into the ocean by rivers is
presented in Fig. <ref>b.
The global, effective river flux of Mn is 0.28.

§.§.§ Sediment source

The largest contribution of Mn to the upper ocean is dust deposition, but over
large shelf and slope regions (e.g. polar oceans) the flux of Mn_diss from the
sediment can be of the same order of magnitude as, or higher than, the dust
deposition flux <cit.>.

The (redox) reactions in the sediment are not explicitly modelled, since
the sediment is not part of our model domain.
Therefore, Mn addition from the sediment is modelled as a prescribed
Since we do not have global maps of Mn_diss sediment–seawater flux, we
parameterise the flux based on existing parameterisations of nitrate and iron.

<cit.> derived an empirical model for calculating the
denitrification rate as a function of the seafloor depth:

\begin{equation}
\zeta_\text{Fsed} = - 0.9543 + 0.7662 \cdot \ln(z_\text{Fsed})
- 0.235 \cdot \ln(z_\text{Fsed})^2 \,,
\end{equation}

where $\zeta_\text{Fsed}$ is the natural logarithm of the denitrification flux
( of carbon) and $z_\text{Fsed}$
is a function of bathymetry.
<cit.> used this model for their sediment source of dissolved Fe in
the Pisces model.
They used a high-resolution bathymetric map to account for the shallow shelves,
and modulated the seafloor depth ($z \rightarrow z_\text{Fsed}$) as follows:

\begin{equation}
z_\text{Fsed} = \min\big(8, \big(\frac{z}{\SI{500}{\metre}}\big)^{-1.5}\big) \,.\\
\end{equation}

We set the maximum iron flux to 1,
following <cit.>.
We assume a porewater ratio $r_\text{Mn:Fe,sed}=0.2$, and follow the same
method as <cit.>.
Our final sediment addition of Mn_diss into bottom water is given

\begin{equation}
\frac{\partial \mathcal{D}}{\partial t}\Big|_\text{sediments} =
%    r_\text{Mn:Fe,sed} \cdot
\frac{\SI{0.2}{\micro\mol\per\square\metre\per\day}}{\Delta z_\text{sed}}
\cdot \min(1,\; 2 \text{e}^{\zeta_\text{Fsed}})\,,
%    \frac{\partial [\chem{Fe}_\text{diss}]}{\partial t}\Big|_\text{sediments}\,,
\label{eqn:mang:sediment}
\end{equation}

where $z_\text{sed}$ is the gridbox thickness of the bottom gridbox (just above
the seafloor) in metre.

With this prescribed source, the Mn flux is limited to a maximum of about
This can be seen in Fig. <ref>c: the
higher-than-70□ regions on the
shelves is especially notable as it falls in the upper (red) part of
the colour scale.
The global sediment flux of Mn is 3.4.

§.§.§ Hydrothermal source

Hydrothermal vent Mn_diss fluxes are modelled proportional to that of ^3He
in <cit.> and <cit.>.
This approach is shown to have worked for iron
The basic equation for the change of Mn_diss from hydrothermal vent influx is

\begin{equation}
\frac{\partial\mathcal{D}}{\partial t}\Big|_\text{hydrothermal} =
r_{\chem{Mn:^3\!He},\text{hydro}} \cdot \frac{\partial [\chem{^3\!He}]}%
{\partial t}\,,
\label{eqn:mang:hydro}
\end{equation}

where $r_\text{\chem{Mn:^3\!He},hydro} = \num{0.1e9}$ is the ratio
between the Mn_diss and ^3He effective inflow
from hydrothermal vents into the model domain.
Recent observational studies found Mn_diss/^3He
concentration ratios of 4.0e6 <cit.> and 3.5e6
<cit.> at observational sites close to hydrothermal outflux regions.
To satisfy this, we assume that a dissolved fraction of 4 % is left when the
hydrothermal plume reaches these two observational sites.

This high “solubility fraction” means that the hydrothermal vents are a large source of
Mn_diss compared to Fe <cit.>.
The integrated Mn flux is 102.
In our simulations we will show that we need to assume such a large flux to
explain the observations.
This choice also relates to the fast modelled removal rate of Mn near the
depth of hydrothermal vents (explained later).
The hydrothermal Mn_diss source is presented in
Fig. <ref>d.

§.§.§ Redox processes

Reduction and oxidation of Mn within the water column is a combination of several
Here we assume that Mn_ox is subject to reduction with a constant
rate, but significantly stimulated by sunlight <cit.>.
This is taken into account in the model by using different $k_\text{red}$ for the
euphotic and aphotic zones of the ocean.
However, there are other processes playing a role that can locally be important.
Those include the microbial enhancement of the rate of oxidation in regions
where Mn_diss supply is high <cit.>, and the
dependence on the O_2 concentration and pH <cit.>.
However, at this stage we decide to not include a dependency on [O_2] to
the model.
Hence, here we choose to model Mn following a pseudo-first-order reaction
where $k_\text{ox}$ and $k_\text{red}$ are pseudo-first-order rate
constants <cit.>:

\begin{align}
\frac{\partial\mathcal{D}}{\partial t}\Big|_{\text{redox}} &=
-k_\text{ox} \mathcal{D} + k_\text{red} \mathcal{X} \\
\frac{\partial\mathcal{X}}{\partial t}\Big|_{\text{redox}} &=
k_\text{ox} \mathcal{D} - k_\text{red} \mathcal{X}  \,,
\label{eqn:redox}
\end{align}

where $\mathcal{X}$ is the particulate oxidised Mn concentration, and

\begin{equation}
k_\text{red} = \begin{cases}
k_\text{red,light} \;\;\;\; \text{in the euphotic zone} \\
k_\text{red,dark} \;\;\;\; \text{elsewhere}    \,,
\end{cases}
\end{equation}

where $k_\text{red,light}$ is the reduction rate in the euphotic zone and
$k_\text{red,dark}$ in the aphotic zone.
The euphotic zone is defined as the depths where the sunlight intensity is at
least 1 %.
The value for $k_\text{red,light}$ is taken from <cit.> who found a
mean dissolution rate of natural Mn oxides of
The $k_\text{red,dark}$ is much smaller and much more uncertain (varying a
factor of fifteen in <cit.>).
As we have both dissolved and particulate Mn in our model, we chose to
fit the $k_\text{ox}/k_\text{red,light}$ to
[Mn_ox]/[Mn_diss] from observations, resulting in a value
within the observational range of <cit.>.
At the time of this study only <cit.> was known to us as reporting
accurate dissolved and particulate Mn for the same samples (at the VERTEX-IV
From this $k_\text{ox}$ was calculated, and as <cit.> sampled the
deep ocean as well, $k_\text{red,dark}$ was derived from that study as well
(Table <ref>).

For one simulation, we will introduce a threshold on the oxidation process of
In that case, oxidation only takes place when [Mn_diss].higher than a certain threshold
value $\mathcal{D}_\text{thr}$; in other words, $k_\text{ox}$ is multiplied by
$H(\mathcal{D} - \mathcal{D}_\text{thr})$.
The Heaviside step function $H(x)$ equals zero where $x\leq 0$ and one where $x>0$.
Based on observations, we have estimated $\mathcal{D}_\text{thr} =
\SI{0.125}{\nano\molar}$, corresponding to the observed deep ocean background
value away from the influence of hydrothermal sources.
The value can be reproduced by first sorting the Mn_diss
Then we cut off the third quartile to remove the high values (above
0.42), and chose a low value for a typical background
concentration (0.13 is the first quartile).
The reason for not choosing simply the minimum value is that values close to that
can be lower than the `background value' because of local removal processes.

§.§.§ Settling and burial

Manganese oxides settle, resulting in a concentration change according to

\begin{equation}
\frac{\partial \mathcal{X}}{\partial t} \Big|_\text{settling}
= - w_\text{ox} \cdot \frac{\partial\mathcal{X}}{\partial z} \,,
\label{eqn:settlingMnox}
\end{equation}

where $w_\text{ox}$ is the settling velocity, set to a constant
1 as long as $\chem{[Mn_\text{ox}]} \leq
\SI{25}{\pico\molar}$.
If $\chem{[Mn_\text{ox}]} > \SI{25}{\pico\molar}$, the settling velocity
is not a constant any more but a function of depth.
Still, in the mixed layer $w_\text{ox} = \SI{1}{\metre\per\day}$, but
if $\chem{[Mn_\text{ox}]} > \SI{25}{\pico\molar}$ the sinking speed
increases linearly such that it reaches a value of 10 at
2.5 depth.
The manganese oxide is buried when arriving at the ocean floor, which
means that it is removed from the model domain.

The rationale for increasing the settling velocity
above a certain threshold is that if Mn_ox is large enough, dense
aggregates form that have faster sinking rates.
This is not unlike increasing velocities of detritus in some models
Notably mineral particles (sand, clay, carbonate) of high density in the
order of 2–3 times that of seawater are responsible for this <cit.>.
Particulate organic carbon may play an important role as well <cit.>.
In this way, [Mn_diss] does not go to zero while still
providing a deep ocean sink.
As long as this critical particulate Mn concentration of 25 is
not reached, aggregation of small Mn_ox particles does not yet
The choice of the critical concentration of 25, also referred
to as the aggregation threshold, is derived from the redox rate constants
in combination with a typical, low value of the observed dissolved Mn
concentration (here chosen as 0.125):
More precisely, settling of Mn_ox follows Eq. (<ref>)
with the settling velocity

\begin{equation}
w_\text{ox}/(\si{\metre\per\day}) = 1 + 9 \cdot
\text{max}\left( 0, \frac{z-\mathit{MLD}}{\SI{2.5}{\kilo\metre}} \right)
\cdot H(\mathcal{X} - \mathcal{X}_\text{thr}) \,,
\label{eqn:settlingMnox_step}
\end{equation}

where $z$ is the depth, MLD is the mixed layer depth, and
$\mathcal{X}_\text{thr}$ is the aggregation threshold set
to 25.

§.§.§ Biological cycle

Manganese is incorporated into phytoplankton during growth.
To this end, we run the Pisces-v2 biogeochemical ocean model
<cit.> together with our manganese model.
The biological processes where Mn is involved are modelled in proportion to the
change in phosphate concentration, thus it is given by:

\begin{equation}
\frac{\partial\mathcal{D}}{\partial t} \Big|_\text{biology} = R_\text{Mn:P}\cdot
\frac{\partial\chem{[PO_4^{3-}]}}{\partial t} \Big|_\text{biology} \,,
\label{eqn:bio}
\end{equation}

where [PO_4^3-] is the phosphate concentration.
The extended Redfield ratio for Mn, $R_\text{Mn:P}$,
is set to 0.36e-3, the value that was determined from data at
the GIPY5 Zero Meridian section <cit.>.
This is a typical value for the manganese to phosphorus ratio in phytoplankton,
though the full range of synchrotron X-ray fluorescence determined ratios,
i.e. those that correspond to intracellular concentrations, is
0.2e-31.5e-3 <cit.>.

There is no growth limitation of phytoplankton by shortage of dissolved Mn,
i.e., in the model, manganese does not affect the biological carbon cycle.

There are four types of plankton in Pisces: nanophytoplankton,
diatoms, microzooplankton and mesozooplankton
(presented in Fig. <ref> as the box “living”).
There are two detrital pools, namely small particles settling with
2 and large particles
settling in our model configuration with 50 (box
Manganese is incorporated in the four living pools according to the
Pisces equations for phosphorus.
Through other processes they become part of detrital material or particulate
organic matter.
These pools are, however, not explicitly followed.
Only the biological sources and sinks of Mn_diss are modelled,
entailing the conversion of Mn_diss to the two phytoplankton
pools, and from all particle pools to Mn_diss.

§.§ Simulations

The Reference simulation was spun up for 600yr to reach a steady state.

From year 100 onwards two sensitivity simulations were forked off to run in
parallel with Reference for another 500.
These simulations are variations of the Reference simulation that uses
the parameters listed in Table <ref>.
Table <ref> lists the simulations and their key parameters.