# Prediction of Tropical Pacific Rain Rates with Over-parameterized Neural
Networks

###### Abstract

The prediction of tropical rain rates from atmospheric profiles poses
significant challenges, mainly due to the heavy-tailed distribution exhibited
by tropical rainfall. This study introduces over-parameterized neural networks
not only to forecast tropical rain rates, but also to explain their heavy-
tailed distribution. The prediction is separately conducted for three rain
types (stratiform, deep convective, and shallow convective) observed by the
Global Precipitation Measurement satellite radar over the West and East
Pacific regions. Atmospheric profiles of humidity, temperature, and zonal and
meridional winds from the MERRA-2 reanalysis are considered as features.
Although over-parameterized neural networks are well-known for their “double
descent phenomenon,” little has been explored about their applicability to
climate data and capability of capturing the tail behavior of data. In our
results, over-parameterized neural networks accurately predict the rain rate
distributions and outperform other machine learning methods. Spatial maps show
that over-parameterized neural networks also successfully describe spatial
patterns of each rain type across the tropical Pacific. In addition, we assess
the feature importance for each over-parameterized neural network to provide
insight into the key factors driving the predictions, with low-level humidity
and temperature variables being the overall most important. These findings
highlight the capability of over-parameterized neural networks in predicting
the distribution of the rain rate and explaining extreme values.

This study aims to introduce the capability of over-parameterized neural
networks, a type of neural network with more parameters than data points, in
predicting the distribution of tropical rain rates from grid-scale
environmental variables and explaining their tail behavior. Rainfall
prediction has been a topic of importance, yet it remains a challenging
problem for its heavy-tailed nature. Over-parameterized neural networks
correctly captured rain rate distributions and the spatial patterns and
heterogeneity of the observed rain rates for multiple rain types, which could
not be achieved by any other previous statistical or machine learning
frameworks. We find that over-parameterized neural networks can play a key
role in general prediction tasks, with potential expanded applicability to
other domains with heavy-tail data distribution.

##  1 Introduction

The tropical Pacific region plays a crucial role in global climate
variability. Rainfall anomalies associated with the El Niño-Southern
Oscillation phenomenon generate atmospheric waves that propagate across the
globe and impact weather in remote regions of the globe. However, climate
models still suffer from biases in this region that affect assessments of
weather and climate risk (Lee et al., 2022; Sobel et al., 2023). Several
attempts have been made to use statistical and machine learning frameworks to
predict rain rates in the tropics. Previous attempts applied machine learning
techniques to the parameterization of convection using cloud-resolving model
or high-resolution model simulations (Brenowitz and Bretherton, 2018; O’Gorman
and Dwyer, 2018; Rasp et al., 2018). These efforts bring computational gains,
but they cannot overcome inherent model deficiencies and may not identify the
underlying physical relationships. Instead, Yang et al. (2019) sought to
capture the relationship between atmospheric features and the rain amount by
applying a generalized linear model (GLM) to observational data. GLMs
naturally provide the relationship through their coefficients. But, they may
struggle to predict the tail of the rain rate distributions accurately due to
their parametric assumption on the density function. Building upon this work,
Wang et al. (2021) improved the predictive performance by implementing random
forest (RF) and neural network (NN) methods. These machine learning methods
allowed for more flexibility of the relationship and demonstrated better
performance than GLM in predicting the rain rate distribution and capturing
its tail behavior. However, the distribution from their outputs still deviated
significantly from the observations. In this work, we experiment with over-
parameterized NNs and find that they enhance the previous results and
effectively capture the tail of the rain rate distribution.

The over-parameterized regime has gained significant attention since the
double descent phenomenon was proposed by Belkin et al. (2019). The double
descent describes the pattern of the test error curve with respect to the
model capacity (Figure 1). According to classical wisdom, the test curve
exhibits a U-shaped pattern before reaching the interpolation threshold, where
perfect interpolation can be attained. Belkin et al. (2019) showed that beyond
the interpolation threshold, the test error decreases again, creating a double
descent pattern of the test curve. Interestingly, the notorious issue of
overfitting, where an excessive number of model parameters compared to the
number of samples harms the test performance, does not show up beyond the
interpolation threshold. Since Belkin et al. (2019), numerous studies have
been conducted to explore and explain this phenomenon theoretically and
experimentally (Geiger et al., 2019; Bartlett et al., 2020; d’Ascoli et al.,
2020; Nakkiran et al., 2020, 2021; Muthukumar et al., 2021; Gamba et al.,
2022; Hastie et al., 2022). However, most of these studies focused on
classification tasks, with little attention given to leveraging the double
descent phenomenon for accurate regression predictions on real datasets. In
addition, to the best of our knowledge, the capability of over-parameterized
NNs in explaining the tail behavior of data has never been investigated. In
this study, we introduce over-parameterized NNs to rainfall predictions and
present numerical findings that demonstrate the efficacy of this approach.

Figure 1: Illustration of the double descent phenomenon. The training curve is
represented by the black line, while the test curve is depicted in blue. The
point on the x-axis where the training error reaches zero is “interpolation
threshold.”

Another crucial aspect of interest in the rainfall dataset is identifying key
features for our model. Even if NNs successfully capture complex relationships
between variables, it may be difficult to disentangle the individual
contributions of each feature. To address this issue, we employed the
permutation importance (PI) approach devised by Breiman (2001). The PI
approach allows us to assess the feature importance by examining the change in
the model output with permuted feature values from the original model output.
By adopting the PI approach, we are able to provide insights into the key
features for our over-parameterized NNs and enhance the interpretability of
our model.

The structure of the remainder of this paper is as follows. We provide a
detailed description of the data used in our study in Section 2. In Section 3,
we recap machine learning methods used in previous studies and introduce the
over-parameterized regime with a suitable training method. Section 4 first
checks if our over-parameterized NNs are properly trained and exhibit the
double descent phenomenon. Then, we compare prediction results of the rain
rate distribution with observations and other machine learning methods.
Lastly, we summarize our findings and discuss future work.

##  2 Data description

We used eight years of June, July and August (JJA) rain rate data from 2015 to
2022 (Figure 2). The rain rate observations were obtained from the Global
Precipitation Measurement (GPM) dual-frequency precipitation radar (DPR)
Version 7 dataset (Iguchi and Meneghini, 2021). JJA data from 2015-2018 were
used for training and JJA data from 2019-2022 were used for testing. The
observational domain was limited to the tropical West Pacific (WP;
$130^{\circ}\text{E}-180^{\circ}\text{E}$,
$15^{\circ}\text{S}-15^{\circ}\text{N}$) and East Pacific (EP;
$180^{\circ}\text{W}-130^{\circ}\text{W}$,
$15^{\circ}\text{S}-15^{\circ}\text{N}$) regions, as illustrated in Figure 2.
These two regions were selected for their distinct convective environments,
with the WP box representing the warm pool region and South Pacific
convergence zone (SPCZ) and the EP box impacted by equatorial upwelling with a
more distinct intertropical convergence zone (ITCZ) north of the equator. The
EP box is also where global climate models (GCMs) commonly experience an
erroneous double ITCZ in boreal summer (e.g., Oueslati and Bellon, 2015). The
DPR data was gridded to a temporal resolution of 3 hours and a spatial
resolution of $0.5^{\circ}$, resulting in 6,000 spatial locations for each
region.

Figure 2: The average GPM DPR rain rate (in mm/day) during JJA 2015-2022 over
the tropical Pacific. The data domains for this study were separated into the
West Pacific (blue) and East Pacific (red) regions.

For the rain types, we adopted a three-type categorization: stratiform, deep
convective, and shallow convective (Funk et al., 2013). This separation allows
us to account for the unique nature of each rain type, with deep convection
extending through the depth of the troposphere due to its stronger updrafts
and shallow convection being confined to echo tops below the $0^{\circ}$C
level because of weaker updrafts and/or a limiting environment (such as dry
air) at upper levels (Schumacher and Funk, 2023). Stratiform rain in the
tropics forms from deep convection and is an important component of mesoscale
convective systems (Houze, 1997). The percent of total rain for each rain type
in each domain is summarized in Table 1. As our focus is on predicting the
rain rate, we only use samples with non-zero rain rates, with numbers of
samples for each region indicated in Table 2. Tables 1 and 2 show that while
the deep convective rain type has the lowest frequencies in both the WP and EP
regions, it contributes about $40\%$ of the total rain. On the contrary, the
shallow convective rain type shows the highest frequencies in both regions,
yet it makes up the smallest portion of the total rain amount (about 10-20%)
in both regions. This emphasizes the distinct nature of each rain type and
demonstrates the importance of separately considering each rain type.

Table 1: Rain amount percentage of each rain type to the total rain amount in the West Pacific and East Pacific regions during JJA. Rain type | WP | EP
---|---|---
Stratiform | 47.2 | 42.1
Deep convective | 42.9 | 37.7
Shallow convective | 9.9 | 20.2
|  |
Table 2: The number of DPR samples with positive rain rate for each rain type and region in the training (JJA 2015-2018) and testing (JJA 2019-2022) datasets. Rain type | WP | EP
---|---|---
Train | Test | Train | Test
Stratiform | 160,170 | 157,682 | 94,034 | 72,898
Deep convective | 113,719 | 115,052 | 67,735 | 52,468
Shallow convective | 207,186 | 208,982 | 185,353 | 170,386
|  |  |  |

As in Yang et al. (2019), we consider humidity, temperature, zonal wind,
meridional wind, latitude, longitude, and latent heat flux as features for the
rain rate prediction. The humidity, temperature, zonal wind, and meridional
wind fields consist of vertical profiles at 40 pressure levels, leading to 163
features in total including latitude, longitude and latent heat flux. The
features used in this study were computed from the MERRA-2 reanalysis
(Rienecker et al., 2011). MERRA-2 data were regridded to have the same
temporal and spatial resolutions as the DPR rain rate dataset.

Rainfall prediction presents a notable challenge due to the extremely thick
tail in the rain rate distribution, indicating a relatively large portion of
extreme values. Table 3 summarizes some statistics of the training dataset for
each rain type and region. The percentiles in Table 3 represent the ranking of
a particular value within a dataset. For instance, the 75th percentile
indicates that 25% of the data values are above it and the 99th percentile
means that only 1% of the data values exceed it. In the case of stratiform and
deep convective rain, the mean exceeds both the median and the 75% percentile
value. This indicates that the data is highly heavy-tailed, as the mean values
are significantly influenced by the presence of extreme values. Although
shallow convective rain may not be as heavy-tailed as the other two rain
types, the mean values lying between the median and 75% percentile values
still points out the presence of a thick tail in the data.

Table 3: Mean, percentiles (median, $75\%$, $90\%$, $99\%$), and maximum DPR rain rate value (in mm/hr) of the training dataset for each rain type in each region for a $0.5^{\circ}$ grid resolution. Region | Rain type | Mean | Median | 75% | 95% | 99% | Max
---|---|---|---|---|---|---|---
WP | Stratiform | 0.42 | 0.10 | 0.31 | 1.86 | 5.42 | 28.74
Deep convective | 0.48 | 0.08 | 0.31 | 2.23 | 6.57 | 60.00
Shallow convective | 0.05 | 0.03 | 0.07 | 0.13 | 0.32 | 2.35
EP | Stratiform | 0.34 | 0.07 | 0.24 | 1.48 | 4.79 | 30.30
Deep convective | 0.47 | 0.10 | 0.34 | 2.15 | 6.18 | 50.38
Shallow convective | 0.08 | 0.04 | 0.10 | 0.19 | 0.48 | 2.97
|  |  |  |  |  |  |

A challenge from dealing with heavy-tailed data arises from a large portion of
extreme samples. These extreme values can distort the overall patterns from
the majority of the data and, hence, decrease the predictive power of models.
In addition, many traditional methods may lack the necessary flexibility to
cover the full range of values from a majority of the data to the extreme
values. As a result, these methods may struggle to estimate and predict the
data, resulting in unsatisfactory overall performance. This limitation can
particularly affect a model’s ability to accurately capture and forecast
extreme events, which is crucial in tasks like rainfall forecasting.

##  3 Methods

In the following, we present four statistical and machine learning frameworks
to investigate the relationship between the rain rates and features.

###  3.1 Generalized linear model

The GLM (McCullagh and Nelder, 1989) is a flexible statistical model that
provides an extension to ordinary linear regression by allowing non-normal
distribution of the response variable. In contrast to ordinary linear
regression, where the mean of the response, $\bm{\mu}$, is directly modeled as
a linear function of the features, the GLM introduces the link function to
establish a nonlinear relationship between the response and the features. The
link function enables the GLM to handle a wider range of data distributions
and provide greater capability in modeling complex relationships.

The GLM consists of three essential components: the random component, the
systematic component, and the link function. The random component deals with
the distributional assumption of the response variable, allowing for non-
normal distributions. The systematic component represents the linear
predictor, which is formed by combining the features with their regression
coefficients. The link function constructs a relationship between the random
and systematic components, i.e.,

$E(Y|\bm{X})=\bm{\mu}=g^{-1}(\bm{X}^{\top}\bm{\beta}),$

where $Y\in\mathbb{R}$ is the response, $\bm{X}\in\mathbb{R}^{p}$ is the
feature with $p$ dimensions, $\bm{\beta}\in\mathbb{R}^{p}$ is the coefficient
of $\bm{X}$, and $g$ is the link function. Depending on the distribution of
the response and the link function, GLMs can take on different names such as
logistic regression, Gamma regression, etc.

Following Yang et al. (2019) and Wang et al. (2021), we implement Gamma
regression for its resemblance to the distribution of rainfall in terms of
skewness. Suppose $Y(\bm{s})$ is a response at location
$\bm{s}\in\mathcal{S}\subset\mathbb{R}^{d}$, where $\mathcal{S}$ is an
observation domain. Then, Gamma regression, with the log-link function
$g(x)=\log(x)$, has the following structure:

$\displaystyle Y(\bm{s})\sim Gamma(\alpha(\bm{s}),\beta(\bm{s})),$
$\displaystyle
E(Y(\bm{s})|\bm{X}(\bm{s}))=\alpha(\bm{s})\cdot\beta(\bm{s})=\mu(\bm{s}),$
$\displaystyle\log(\mu(\bm{s}))=\eta_{0}+\eta_{1}X_{1}(\bm{s})+\cdots\eta_{p}X_{p}(\bm{s}),$

where $\alpha(s)>0$ and $\beta(s)>0$ are parameters from the Gamma
distribution where a probability density function is given as

$f(y;\alpha,\beta)=\displaystyle\frac{1}{\beta^{\alpha}\Gamma(\alpha)}y^{\alpha-1}e^{-\frac{y}{\beta}},~{}~{}0<y<\infty.$

Here, $\Gamma$ is a Gamma function, $\bm{X}(\bm{s})$ is the feature vector at
location $\bm{s}$, and $\eta_{i}\text{'s}\ (i=0,\ldots,p)$ are coefficients
for each feature with the total number of $p$. Parameter estimation in Gamma
regression can be done by maximum likelihood estimation, which is provided by
R with glm function.

Despite some advantages of Gamma regression, such as its ability to handle
skewed data, its limited parametric distribution we assume for the error may
not adequately capture the characteristics of heavy-tailed data such as
rainfall (Yang et al., 2019; Wang et al., 2021). Thus, adopting machine
learning techniques, which do not impose such strong restrictions as GLMs,
could potentially enhance the prediction performance.

###  3.2 Random forest

RF has gained attention for its powerful prediction capability since it first
appeared in Breiman (2001). RF is an ensemble method that aggregates multiple
decision trees to achieve accurate and robust predictions. A decision tree is
a non-parametric machine learning method that has a hierarchical tree
structure. It consists of three key components: internal nodes, branches, and
leaf nodes. Each internal node represents a condition on one of the features.
The branches indicate split data based on the outcomes of the condition at
each internal node. The leaf nodes denote the final outcomes from split data
at corresponding branches.

Ensemble methods summarize the results from multiple statistical or machine
learning methods into one final outcome. In an RF model, each decision tree
independently learns from random samples of training data with a random subset
of features at each split. This random sampling and subsetting process help to
bring diversity and reduce the risk of overfitting by preventing the trees
from becoming too specialized to the training data or correlated to one
another. Aggregating the predictions from multiple less correlated trees
results in outcomes with lower variance and improved accuracy.

Although RF may lack interpretability compared to GLM, its performance in
fitting and predicting rainfall data exceeded GLM and matched NN (Wang et al.,
2021). In our study, we train a random forest model in R with randomForest
function supported by randomForest package. The number of trees was set to 100
and the number of features selected at each split was set to the largest
number below $\sqrt{\text{the total number of features}}$, which is 12. For
the other configurations, we used the default values provided with
randomForest function.

###  3.3 Neural networks

####  3.3.1 Multi-layer perceptron NN

As the size of datasets grow and high-performing computation becomes more
available, NNs have been one of the most popular frameworks in the statistical
and machine learning communities. Inspired by the structure of the human
brain, NNs are computational models which are defined by hidden units and
connections between the hidden units. Hidden units and connections correspond
to the neurons and the synapses in human brain, respectively. A typical
example would be multi-layer perceptron (MLP), which connects all the units
from the input to the output (LeCun et al., 2015). Besides MLP, convolutional
NN (Simonyan and Zisserman, 2014), recurrent NN (Hochreiter and Schmidhuber,
1997), and transformer (Vaswani et al., 2017) are frequently used structures.
In this study, MLP is used to compare standard NN performance with GLM and RF,
as our primary goal lies in introducing over-parameterized NNs.

MLP is usually established with an input layer followed by sequential hidden
layers and an output layer at the end. The number of layers is often called
the depth of the NN and the number of units in each layer is called the width
of the NN. To take into account nonlinear structure, nonlinear activation
functions are applied after each hidden layer. Some popular nonlinear
activation functions are Tangent hyperbolic (Tanh), Rectified Linear Unit
(ReLU), Gaussian Error Linear Unit (GELU), etc. The choice of the nonlinear
activation function often plays a key role in the model performance and
optimization. In this study, we chose either the Tanh or ReLU function that
best explains each rain type in terms of test errors. At the output layer, an
exponential function is applied to guarantee positive values of the outputs.
In short, MLP with $L$ layers can be put down in the following equations:

Input layer: $\displaystyle\bm{Z}^{(1)}=\begin{bmatrix}Z_{1}\\\ Z_{2}\\\
\vdots\\\ Z_{p}\end{bmatrix}$ Hidden layer 1:
$\displaystyle\bm{H}^{(1)}=\bm{W}^{(1)}\bm{Z}^{(1)}+\bm{b}^{(1)}$
$\displaystyle\bm{Z}^{(2)}=\sigma^{(1)}(\bm{H}^{(1)})$ Hidden layer 2:
$\displaystyle\bm{H}^{(2)}=\bm{W}^{(2)}\bm{Z}^{(2)}+\bm{b}^{(2)}$
$\displaystyle\bm{Z}^{(3)}=\sigma^{(2)}(\bm{H}^{(2)})$ $\displaystyle\vdots$
Output layer:
$\displaystyle\bm{H}^{(L-1)}=\bm{W}^{(L-1)}\bm{Z}^{(L-1)}+\bm{b}^{(L-1)}$
$\displaystyle Y=\sigma^{(L-1)}(\bm{H}^{(L-1)}),$

where $\bm{Z}^{(l)}$ is the input vectors to $l$-th layer, $\bm{H}^{(l)}$ is
the hidden outputs from $l$-th linear layer and $Y$ is the final output.
$\bm{W}^{(l)}$’s and $\bm{b}^{(l)}$’s are the weights and biases of the NNs
and $\sigma^{(l)}$ is the activation function at $l$-th hidden layer.

####  3.3.2 Over-parameterized regime

An over-parameterized regime refers to when the number of parameters in a
model exceeds the training sample size. Over-parameterized models have gained
significant attention in recent years (Cao et al., 2022; Hastie et al., 2022;
Liu et al., 2022; Xu et al., 2022; Gao et al., 2023), particularly with the
rise of deep learning, as many state-of-the-art deep learning methods fall
into this category. From classical wisdom, an excess number of parameters
compared to the number of data points usually leads to overfitting. However,
it has been observed that over-parameterized NN models do not suffer from
traditional overfitting and result in decent generalization.

Belkin et al. (2019) experimentally showed how classical understanding of
overfitting may have to change in an over-parameterized regime. Their study
demonstrated that over-parameterized models not only perfectly interpolate the
training data, but also attain lower test error than what can be achieved in
under-parameterized regimes. According to the classical understanding, the
test error curve first decreases until an unknown “sweet spot” where the
minimum test error is obtained, then it starts to increase again. The maximum
test error occurs around the interpolation threshold, where the number of
parameters matches the sample size. Beyond the interpolation threshold, also
called the over-parameterized regime, the test error starts to decrease again.
This phenomenon is often called double descent.

However, it becomes more challenging to train over-parameterized models as we
need to explore a larger parameter space. Belkin et al. (2019) suggested a
sequential training procedure to gradually increase the model size, with the
final weight of smaller models serving as the pre-trained initialization for
the larger models. We followed a similar sequential training procedure and
considered NNs with different depths ($L=4,5,6$) and widths
($W=12,24,\ldots,600$). For example, when training NNs with $L=4$, the learned
weights of a NN with $L=4$ and $W=12$ were used as the initialization for a NN
with $L=4$ and $W=24$. The remaining weights were initialized by default in
PyTorch, which in our case was the Kaiming uniform initialization (He et al.,
2015). For deeper NNs ($L=5,6$), the same width was maintained, but NNs with
one less layer were used for initialization. Our final models were NNs with
$L=6$, but NNs with $L=5$ were considered for initialization of NNs with
$L=6$, as we found this process stabilized the training procedure. This pre-
training strategy allowed us to optimize over a larger parameter space while
benefiting from the knowledge gained by smaller models in the sequential
progression.

Every experiment with NNs was conducted by PyTorch with the Adam optimizer
(Kingma and Ba, 2014). The total number of epochs for every NN with a
different size was set to be 2000 and the learning rate decayed by 0.9 every
100 epochs. To quantify the discrepancy between predicted and target values,
we used the mean squared error (MSE) as the loss function. It is important to
note that no explicit or implicit regularization techniques, such as dropout
or imposing a penalty on weights, were utilized. More details on the training
configurations for each region and rain type are listed in Table 4.

Table 4: Training configurations for NNs with 6 layers. WP and EP stand for the West Pacific and East Pacific region, respectively. Region | Rain type | Activation | Initial learning rate
---|---|---|---
WP | Stratiform | Tanh | $5\times 10^{-5}$
Deep convective | Tanh | $9\times 10^{-5}$
Shallow convective | Tanh | $10\times 10^{-5}$
EP | Stratiform | Tanh | $10\times 10^{-5}$
Deep convective | ReLU | $1\times 10^{-5}$
Shallow convective | Tanh | $10\times 10^{-5}$
|  |  |

###  3.4 Model interpretation

One of the popular approaches to explain variable importance in machine
learning is permutation importance (Breiman, 2001; Fisher et al., 2019; Sood
and Craven, 2022; Cheung et al., 2022; Ramos-Valle et al., 2023; Sekiyama et
al., 2023). Permutation importance quantifies the importance of individual
features by assessing the impact on test errors when their values are randomly
shuffled while keeping the other features unchanged. To calculate the
permutation importance for a specific feature (e.g., temperature at a
particular height), we first shuffle values of the target feature while
leaving the remaining features in place. Then, we obtain model predictions
from this permuted data and compute the increase in the loss values. A larger
increase in loss values indicates that the target feature plays a more
significant role in the model, as breaking the relationship between this
target feature and the model outcome leads to a greater deterioration in model
performance.

The original idea, proposed by Breiman (2001), aimed to explore local variable
importance, i.e., the variable importance for an individual sample. However,
our primary objective is to achieve transparency for the entire model, hence
we employ permutation importance to derive global variable importance (König
et al., 2021; Sood and Craven, 2022). The only distinction between calculating
local and global variable importance lies in whether each permutation is kept
distinct or identical across samples for a given target feature. While there
are other popular interpretation methods such as SHapley Additive exPlanations
(SHAP, Lundberg and Lee (2017)) and Local Interpretable Model-agnostic
Explanations (LIME, Ribeiro et al. (2016)), they mainly provide local variable
importance. Since our focus is on global variable importance, permutation
importance was considered the most suitable choice for this study.

To summarize our variable importance findings, we employed the following
process. For each feature, we conducted 100 random permutations, ensuring the
robustness of our variable importance investigations. In each permutation, we
identified the top 10 most significant features and aggregated these
selections across all permutations. Finally, we determined key features for
our model by reporting the top 10 most frequently occurring ones within this
combined set of features.

##  4 Results

Only data points with positive rain rates were used for the results presented
below. We compare four statistical and machine learning methods: GLM, RF, the
under-parameterized NN, and the over-parameterized NN. We chose the under-
parameterized NN model that yielded the best test error before peaks in the
test curve and the over-parameterized NN model with the largest number of
parameters. The best models were separately obtained for the WP and EP regions
and three different rain types: stratiform, deep convective, and shallow
convective.

###  4.1 Double descent

We first examine whether over-parameterized NNs are well trained. To properly
train over-parameterized NNs, the model size should be increased beyond the
interpolation threshold, where we can achieve zero training error and observe
the double descent curve (Belkin et al., 2019, Figure 1). Then, over-
parameterized models at the end of the double descent curve are used to see
how these models perform with heavy-tailed rainfall data.

Figure 3 illustrates the training and testing curves obtained from NNs with
varying depths for different regions and rain types. As the model size
increases, we are able to achieve nearly zero training error across all rain
types. However, due to the heavy-tailed nature of the rainfall dataset, larger
models were necessary to attain sufficiently low training error beyond the
typical interpolation thresholds (i.e., where the sample size is equal to the
number of model parameters as indicated by the vertical dashed line in each
panel). In particular, NNs for deep convective rain need larger models to
achieve nearly zero training error, as it is the most heavy-tailed rain type.
Nakkiran et al. (2021) also pointed out that other factors such as data
distribution and training procedure may have influence on the model
complexity, which is directly related to the double descent phenomenon.

Notably, all the test curves in Figure 3 exhibit the double descent
phenomenon. The test errors after the peaks decrease and reach minimum or
near-minimum points, although the peaks in the test curves occur at different
positions. This observation holds true for all regions and rain types examined
in the study, indicating that the over-parameterized NNs have been
legitimately trained. A minor point to discuss is the position of the peaks in
the test errors. In theory, with simple frameworks such as linear regression,
the peaks are expected to occur at the interpolation threshold. However, a
shift in the peak position has been observed in the literature and various
reasons such as the data distribution, noises in responses, or regularization
may cause the shift (d’Ascoli et al., 2020; Nakkiran et al., 2021; Gamba et
al., 2022).

Figure 3: Training and test curves for the West Pacific (left two columns) and
East Pacific (right two columns). The training results are the left panel and
the test results are the right panel in each subfigure. The ratio of the
number of parameters to the number of samples is on the x-axis, and the RMSE
values are on the y-axis. The theoretical interpolation threshold is marked as
a blue dashed line. Black, red, and blue curves represent NNs with 4, 5, and 6
layers, respectively.

###  4.2 Rain rate distributions

Figure 4 illustrates percentile plots of estimated (Train) and predicted
(Test) rain rate values from the GLM, RF, and NN methods over the WP and EP.
For both regions and all rain types, the distributions of estimated values
from the over-parameterized NNs (red dashed line) successfully align with the
distribution of the DPR observations (solid black line). This result is
expected given that over-parameterized NNs achieved nearly zero training error
in Figure 3. The other methods could not fully recover the training data and
overestimate occurrence at low rain rates and underestimate occurrence at
large rain rates. This result is consistent with Yang et al. (2019) and Wang
et al. (2021).

In the prediction results, over-parameterized NNs outperform the other methods
and accurately capture the general shape of the rain rate distributions of the
test dataset (Figure 4). While the test lines for the over-parameterized NNs
no longer lie directly on the observed DPR rain rates, they still hew closely
to each other at both low to moderate percentiles (0 to 90%) and tail values
(90 to 99.9%). There are generally better rain rate predictions over the EP,
similar to the GLM performance in Yang et al. (2019). The deep convective
prediction has especially high fidelity in the EP, whereas deep convective
rain rates are all overpredicted in the WP.

Figure 4 also shows that all the other methods exhibit markedly different
patterns in their prediction distributions compared to the true rain rate
distribution. Both GLM and RF tend to overpredict small rain rate values in
the test data while underpredicting large values. This behavior is often seen
in GCMs (e.g., Dai, 2006; Fiedler et al., 2020), and indicates that these
methods produce a narrower range of rain rate predictions compared to the
wider distributions in over-parameterized NNs and the actual observations.
Under-parameterized NNs initially show an underprediction at low rain rates,
but then they bounce back to the observed rain rate distribution at higher
values. However, under-parameterized NNs eventually end up generating too
large of rain rates compared to the true values. Overall, these findings
verify the superiority of the over-parameterized NNs in accurately predicting
rain rate distributions.

Figure 4: Rain rate percentile plots (in mm/hr) for $0.5^{\circ}$ grids in the
West Pacific and East Pacific regions for each rain type. Obs, Under NN, and
Over NN stand for DPR observations, under-parameterized NN, and over-
parameterized NN, respectively. Values on the x-axis indicate the 90%, 99% and
99.9% percentiles from the training and test datasets for each domain and rain
type.

###  4.3 Spatial maps

We next present the geographical patterns produced by each method. The rain
rate prediction results are averaged over multiple years (2019-2022) to
generate a single spatial map for each method. Figures 5 and 6 display
averaged spatial maps over the WP and EP regions, respectively. While the
observations indicate strong spatial heterogeneity over the domains for each
rain type, RF and GLM produce very smooth prediction maps. Under-parameterized
NNs produce non-smooth spatial patterns, but the high rain rates are too
clustered and large. Over-parameterized NNs, on the other hand, demonstrate
their remarkable ability to capture the spatial heterogeneity of the DPR rain
rate observations. Their pebble-like quality is much more representative of
true weather events, compared to the under-parameterized NN boulders and
ruffled sand of the RF and GLM maps.

In terms of overall spatial patterns, all the methods capture the widespread
rainfall of the WP warm pool (Figure 5) and the ITCZ in the northern part of
the EP domain (Figure 6). However, the over-parameterized NN produces higher
fidelity in the patterns. For example, while most of the training was done
over ocean grid points, the over-parameterized NN provides more distinct
predictions over Papua New Guinea in the southwest corner of the WP domain and
is able to capture orographic features to a higher degree than the other
methods, especially in the shallow convective field (Figure 5c). This result
is promising for the future use of over-parameterized NNs for rainfall
prediction over continents. In the EP, the over-parameterized NN produces a
more distinct ITCZ in the stratiform and deep convective fields (Figures 6a
and b), while also reproducing the pattern in shallow convective rainfall at
the equator (Figure 6c) despite the fact that sea surface temperatures are not
used as a predictor.

Figure 5: Maps of the predicted mean $0.5^{\circ}$ JJA rain rates (in mm/hr)
from 2019-2022 over the West Pacific region for (a) stratiform, (b) deep
convective, and (c) shallow convective rain types using results from the over-
parameterized NN, under-parameterized NN, RF, and GLM. Figure 6: As in Figure
5, but over the East Pacific region.

###  4.4 Interpreting over-parameterized neural networks

In this section, we provide the feature importance for over-parameterized NNs
with permutation importance. Table 5 presents the 10 top features for the test
output, categorized into five groups: humidity, temperature, zonal wind,
meridional wind, and others. Each category comprises values at 40 pressure
levels except for the ‘Others’ category, which contains fields reported at
only one level (i.e., longitude, latitude, and latent heat flux). The
individual feature contributions to the model output might not hold
significant meaning, since permuting one variable while keeping the other
variables the same may lead to unrealistic environmental conditions. Instead,
our primary focus lies in identifying which category plays a significant role
in the model output. Future work would be to analyze the individual and
combined feature contributions in more detail.

Higher numbers in Table 5 represent a larger importance. For example, six
pressure levels from the humidity category at or below 750 hPa are identified
as the most significant for predicting shallow convective rain rates in the
EP, while the other four most important features are temperature values at low
(925 and 875 hPa) and mid-to-upper (400 and 350 hPa) levels. While we used all
available height levels (including the stratosphere) from MERRA-2 in our
training and testing, we only report on the importance of pressure levels
within the troposphere (i.e., up to 100 hPa), as almost all weather phenomena
occur within this atmospheric layer.

Table 5: Categories of top 10 important features from test data for the West Pacific and East Pacific regions. The numbers in parentheses represent pressure levels (in hPa) of selected features.  Region | Rain type | Humidity | Temperature | Zonal wind | Meridional wind | Others |
---|---|---|---|---|---|---|---
|  |  |  |  |  |  |
| Stratiform | 5 | 4 | 0 | 1 | 0 |
|  | (875 800 775 725 500) | (1,000 925 800 650) |  | (300) |  |
WP | Deep convective | 4 | 4 | 1 | 0 | 1 |
|  | (875 700 550 450) | (900 800 775 700) | (550) |  | (Latitude) |
| Shallow convective | 5 | 3 | 0 | 2 |  |
|  | (1,000 975 950 400 300) | (1,000 975 300) |  | (975 250) |  |
|  |  |  |  |  |  |
| Stratiform | 2 | 5 | 2 | 0 | 1 |
|  | (1,000 950) | (975 950 900 875 850) | (825 450) |  | (Latitude) |
EP | Deep convective | 4 | 4 | 2 | 0 | 0 |
|  | (725 700 600 350) | (925 850 400 175) | (125 100) |  |  |
| Shallow convective | 6 | 4 | 0 | 0 | 0 |
|  | (1,000 975 950 925 775 750) | (925 875 400 350) |  |  |  |
Average | 4.33 | 4.00 | 0.83 | 0.50 | 0.33 |
|  |  |  |  |  |  |

Table 5 shows that the average number of variables in the humidity category
that are selected as key features is 4.33, indicating that 4.33 variables from
the humidity category were consistently chosen as the 10 most influential
features in rain rate prediction across the tropical Pacific. This result is
consistent with our common understanding that humidity has a direct
relationship with rain rate (e.g., Bretherton et al., 2004), making it a
critical factor in predicting rainfall. Furthermore, for most of the rain
types, humidity contains the most variables among key features, which is
consistent with the results of Ahmed and Schumacher (2015) that indicated a
strong relationship between column humidity and shallow convective, deep
convective, and stratiform rain rates across the tropical oceans. The humidity
category is closely followed by temperature in importance, which aligns well
with the results from Yang et al. (2019). We note that most of the pressure
levels highlighted for humidity and temperature are at or below 700 hPa.

The most influential categories in explaining the rainfall prediction results
over the tropical Pacific in Table 5 after humidity and temperature are, in
their respective order, zonal wind, meridional wind, and others. Even though
wind (or wind shear) is not typically included as a parameter in GCM
convective parameterizations, the small but perceptible presence of zonal and
meridional wind in Table 5 suggests that the inclusion of winds has the
potential to improve convective and stratiform rain predictions. More precise
quantification of these relationships are left for future studies, however,
differences in the category rankings between the WP and EP and between rain
types are consistent with our general understanding of large-scale
environmental factors that affect rain production in the tropics.

##  5 Discussion

Capturing the observed climate variability in the tropical Pacific is one of
the great challenges currently faced by climate modeling (Lee et al., 2022;
Sobel et al., 2023). The rainfall response to sea surface temperature
anomalies in this region plays a key role in driving this climate variability.
In this study, we applied statistical and machine learning techniques to model
the relationship between the large-scale environment and rainfall. With
growing interest in the application of machine learning to climate modeling
(Brenowitz and Bretherton, 2018; O’Gorman and Dwyer, 2018; Rasp et al., 2018),
our results can help guide the development of new machine learning-derived
parameterizations of rainfall.

We found that properly trained over-parameterized NNs correctly predicted the
rain rate distributions over the tropical Pacific for multiple rain types
(stratiform, deep convective, and shallow convective), including their tail
behavior. Over-parameterized NNs also outperformed the other methods such as
GLM and RF in both predicting the rain rate distribution and capturing spatial
patterns. The permutation importance was implemented to address the feature
importance in the NN model outcomes, producing consistent results with those
obtained from the GLM framework in Yang et al. (2019), who also found humidity
and temperature to be the most important environmental variables. To the best
of our knowledge, the benefits of over-parameterized models have been neither
experimentally nor theoretically verified for heavy-tailed datasets. Given our
successful results, the applicability of over-parameterized NNs can be
expanded to various other real-world datasets.

It is notable that the training and test curves in Figure 3 show slightly
different behaviors compared to previous studies. The interpolation thresholds
have shifted to the right due to the heavy-tailed nature of the dataset, as we
need larger model complexity to fully explain samples with large rain rate
values. The locations of the peaks in the test curves were often placed much
earlier than the interpolation thresholds. While it is partly known that
factors like data distribution and noise in responses can influence the
locations of the peaks, the majority of reasons and their impacts remain
unclear and unexplored. One workaround with this issue is to consider
$\ell_{2}$-regularization to the weight of NNs (Nakkiran et al., 2020). They
found that with the optimally chosen tuning parameter, the peak can be
completely eliminated, and even non-optimal regularization can alleviate the
peak. Additionally, optimal early stopping based on test errors may also
contribute to removing the peak, though this only happened for certain cases
(Nakkiran et al., 2021). Our study does not delve further into this matter, as
our primary goal is to introduce over-parameterized NNs to rainfall prediction
and examine the prediction performance of the method. Our findings indicate
that no matter where the peaks are observed, we can attain minimum or near-
minimum test errors in the over-parameterized regime and the obtained models
correctly predict the rain rate distribution. We hope that future research,
both in theory and applications, will shed light on the different patterns of
the double descent curve that we observed with heavy-tailed data.

In addition, the final models for each rain type and region have the potential
for further improvement through fine-tuning and optimal selection of the
hyperparameters. For instance, as observed in Figure 3, deeper neural networks
consistently lead to better prediction performance even when the number of
parameters is similar. The choice of activation functions and the optimizer
also played a crucial role in determining the model performance. However, the
process of searching for optimal hyperparameters needs to be done for each
dataset individually, and it can be tedious work. Furthermore, the
improvements achieved through hyperparameter tuning may not yield significant
differences in the final results. We followed Nakkiran et al. (2021) to choose
the Adam optimizer over stochastic gradient descent and used an MLP structure
to claim that we can obtain desirable results even with simple models. While
further fine-tuning and customization of the model structure could potentially
lead to improvements, the current method allows us to obtain promising results
without the exhaustive search for hyperparameters. This helps maintain a
balance between model performance and computational efficiency, making the
approach more practical and applicable to real datasets.

###### Acknowledgements.

Aaron Funk processed the gridded DPR and MERRA-2 data used in the analysis.
Mikyoung Jun acknowledges support by NSF IIS-2123247 and DMS-2105847. Courtney
Schumacher acknowledges support by NASA 80NSSC22K0617. Funding for this
project was also provided by the TAMU College of Arts & Sciences Seed Funds.
The GPM DPR and MERRA-2 data are publicly available from the NASA GES DISC
(https://disc.gsfc.nasa.gov/). The output from the statistical and machine
learning analysis of the data are available upon reasonable request from the
authors.

## References

  * Ahmed and Schumacher (2015) Ahmed, F., and C. Schumacher, 2015: Convective and stratiform components of the precipitation-moisture relationship. Geophysical Research Letters, 42 (23), 10–453.
  * Bartlett et al. (2020) Bartlett, P. L., P. M. Long, G. Lugosi, and A. Tsigler, 2020: Benign overfitting in linear regression. Proceedings of the National Academy of Sciences, 117 (48), 30 063–30 070.
  * Belkin et al. (2019) Belkin, M., D. Hsu, S. Ma, and S. Mandal, 2019: Reconciling modern machine-learning practice and the classical bias–variance trade-off. Proceedings of the National Academy of Sciences, 116 (32), 15 849–15 854.
  * Breiman (2001) Breiman, L., 2001: Random forests. Machine learning, 45, 5–32.
  * Brenowitz and Bretherton (2018) Brenowitz, N. D., and C. S. Bretherton, 2018: Prognostic validation of a neural network unified physics parameterization. Geophysical Research Letters, 45 (12), 6289–6298.
  * Bretherton et al. (2004) Bretherton, C. S., M. E. Peters, and L. E. Back, 2004: Relationships between water vapor path and precipitation over the tropical oceans. Journal of climate, 17 (7), 1517–1528.
  * Cao et al. (2022) Cao, Y., Z. Chen, M. Belkin, and Q. Gu, 2022: Benign overfitting in two-layer convolutional neural networks. Advances in neural information processing systems, 35, 25 237–25 250.
  * Cheung et al. (2022) Cheung, H. M., C.-H. Ho, and M. Chang, 2022: Hybrid neural network models for postprocessing medium-range forecasts of tropical cyclone tracks over the western north pacific. Artificial Intelligence for the Earth Systems, 1 (4), e210 003.
  * Dai (2006) Dai, A., 2006: Precipitation characteristics in eighteen coupled climate models. Journal of climate, 19 (18), 4605–4630.
  * d’Ascoli et al. (2020) d’Ascoli, S., M. Refinetti, G. Biroli, and F. Krzakala, 2020: Double trouble in double descent: Bias and variance (s) in the lazy regime. International Conference on Machine Learning, PMLR, 2280–2290.
  * Fiedler et al. (2020) Fiedler, S., and Coauthors, 2020: Simulated tropical precipitation assessed across three major phases of the coupled model intercomparison project (cmip). Monthly Weather Review, 148 (9), 3653–3680.
  * Fisher et al. (2019) Fisher, A., C. Rudin, and F. Dominici, 2019: All models are wrong, but many are useful: Learning a variable’s importance by studying an entire class of prediction models simultaneously. J. Mach. Learn. Res., 20 (177), 1–81.
  * Funk et al. (2013) Funk, A., C. Schumacher, and J. Awaka, 2013: Analysis of rain classifications over the tropics by version 7 of the trmm pr 2a23 algorithm. Journal of the Meteorological Society of Japan. Ser. II, 91 (3), 257–272.
  * Gamba et al. (2022) Gamba, M., E. Englesson, M. Björkman, and H. Azizpour, 2022: Deep double descent via smooth interpolation. Transactions on Machine Learning Research.
  * Gao et al. (2023) Gao, Y., S. Mahadevan, and Z. Song, 2023: An over-parameterized exponential regression. arXiv preprint arXiv:2303.16504.
  * Geiger et al. (2019) Geiger, M., S. Spigler, S. d’Ascoli, L. Sagun, M. Baity-Jesi, G. Biroli, and M. Wyart, 2019: Jamming transition as a paradigm to understand the loss landscape of deep neural networks. Physical Review E, 100 (1), 012 115.
  * Hastie et al. (2022) Hastie, T., A. Montanari, S. Rosset, and R. J. Tibshirani, 2022: Surprises in high-dimensional ridgeless least squares interpolation. Annals of statistics, 50 (2), 949.
  * He et al. (2015) He, K., X. Zhang, S. Ren, and J. Sun, 2015: Delving deep into rectifiers: Surpassing human-level performance on imagenet classification. Proceedings of the IEEE international conference on computer vision, 1026–1034.
  * Hochreiter and Schmidhuber (1997) Hochreiter, S., and J. Schmidhuber, 1997: Long short-term memory. Neural computation, 9 (8), 1735–1780.
  * Houze (1997) Houze, R. A. J., 1997: Stratiform precipitation in regions of convection: A meteorological paradox? Bulletin of the American Meteorological Society, 78 (10), 2179–2196.
  * Iguchi and Meneghini (2021) Iguchi, T., and R. Meneghini, 2021: Gpm dpr precipitation profile l2a 1.5 hours 5 km v07. Greenbelt, MD, Goddard Earth Sciences Data and Information Services Center (GES DISC), Accessed, 10.5067/GPM/DPR/GPM/2A/07.
  * Kingma and Ba (2014) Kingma, D. P., and J. Ba, 2014: Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980.
  * König et al. (2021) König, G., C. Molnar, B. Bischl, and M. Grosse-Wentrup, 2021: Relative feature importance. 2020 25th International Conference on Pattern Recognition (ICPR), IEEE, 9318–9325.
  * LeCun et al. (2015) LeCun, Y., Y. Bengio, and G. Hinton, 2015: Deep learning. nature, 521 (7553), 436–444.
  * Lee et al. (2022) Lee, S., M. L’Heureux, A. T. Wittenberg, R. Seager, P. A. O’Gorman, and N. C. Johnson, 2022: On the future zonal contrasts of equatorial pacific climate: Perspectives from observations, simulations, and theories. npj Clim Atmos Sci, 5 (82).
  * Liu et al. (2022) Liu, C., L. Zhu, and M. Belkin, 2022: Loss landscapes and optimization in over-parameterized non-linear systems and neural networks. Applied and Computational Harmonic Analysis, 59, 85–116.
  * Lundberg and Lee (2017) Lundberg, S. M., and S.-I. Lee, 2017: A unified approach to interpreting model predictions. Advances in neural information processing systems, 30.
  * McCullagh and Nelder (1989) McCullagh, P., and J. A. Nelder, 1989: Generalized linear models. 2nd ed., Chapman and Hall.
  * Muthukumar et al. (2021) Muthukumar, V., A. Narang, V. Subramanian, M. Belkin, D. Hsu, and A. Sahai, 2021: Classification vs regression in overparameterized regimes: Does the loss function matter? The Journal of Machine Learning Research, 22 (1), 10 104–10 172.
  * Nakkiran et al. (2021) Nakkiran, P., G. Kaplun, Y. Bansal, T. Yang, B. Barak, and I. Sutskever, 2021: Deep double descent: Where bigger models and more data hurt. Journal of Statistical Mechanics: Theory and Experiment, 2021 (12), 124 003.
  * Nakkiran et al. (2020) Nakkiran, P., P. Venkat, S. Kakade, and T. Ma, 2020: Optimal regularization can mitigate double descent. arXiv preprint arXiv:2003.01897.
  * O’Gorman and Dwyer (2018) O’Gorman, P. A., and J. G. Dwyer, 2018: Using machine learning to parameterize moist convection: Potential for modeling of climate, climate change, and extreme events. Journal of Advances in Modeling Earth Systems, 10 (10), 2548–2563.
  * Oueslati and Bellon (2015) Oueslati, B., and G. Bellon, 2015: The double itcz bias in cmip5 models: Interaction between sst, large-scale circulation and precipitation. Climate dynamics, 44, 585–607.
  * Ramos-Valle et al. (2023) Ramos-Valle, A. N., J. Alland, and A. Bukvic, 2023: Using machine learning to understand relocation drivers of urban coastal populations in response to flooding. Artificial Intelligence for the Earth Systems, 2 (2), 220 054.
  * Rasp et al. (2018) Rasp, S., M. S. Pritchard, and P. Gentine, 2018: Deep learning to represent subgrid processes in climate models. Proceedings of the National Academy of Sciences, 115 (39), 9684–9689.
  * Ribeiro et al. (2016) Ribeiro, M. T., S. Singh, and C. Guestrin, 2016: ” why should i trust you?” explaining the predictions of any classifier. Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining, 1135–1144.
  * Rienecker et al. (2011) Rienecker, M. M., and Coauthors, 2011: Merra: Nasa’s modern-era retrospective analysis for research and applications. Journal of climate, 24 (14), 3624–3648.
  * Schumacher and Funk (2023) Schumacher, C., and A. Funk, 2023: Assessing convective-stratiform precipitation regimes in the tropics and extratropics with the gpm satellite radar. Geophysical Research Letters, 50 (14), e2023GL102 786.
  * Sekiyama et al. (2023) Sekiyama, T. T., S. Hayashi, R. Kaneko, and K.-i. Fukui, 2023: Surrogate downscaling of mesoscale wind fields using ensemble super-resolution convolutional neural networks. Artificial Intelligence for the Earth Systems, 1–35.
  * Simonyan and Zisserman (2014) Simonyan, K., and A. Zisserman, 2014: Very deep convolutional networks for large-scale image recognition. arXiv preprint arXiv:1409.1556.
  * Sobel et al. (2023) Sobel, A. H., and Coauthors, 2023: Near-term tropical cyclone risk and coupled earth system model biases. Proceedings of the National Academy of Sciences, 120 (33), e2209631 120.
  * Sood and Craven (2022) Sood, A., and M. Craven, 2022: Feature importance explanations for temporal black-box models. Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 36, 8351–8360.
  * Vaswani et al. (2017) Vaswani, A., N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, 2017: Attention is all you need. Advances in neural information processing systems, 30.
  * Wang et al. (2021) Wang, J., R. K. Wong, M. Jun, C. Schumacher, R. Saravanan, and C. Sun, 2021: Statistical and machine learning methods applied to the prediction of different tropical rainfall types. Environmental research communications, 3 (11), 111 001.
  * Xu et al. (2022) Xu, W., R. T. Chen, X. Li, and D. Duvenaud, 2022: Infinitely deep bayesian neural networks with stochastic differential equations. International Conference on Artificial Intelligence and Statistics, PMLR, 721–738.
  * Yang et al. (2019) Yang, J., M. Jun, C. Schumacher, and R. Saravanan, 2019: Predictive statistical representations of observed and simulated rainfall using generalized linear models. Journal of climate, 32 (11), 3409–3427.


find $\Delta S_{8}=0.73\sigma$ when comparing the two HMCode2016 analyses and
$\Delta S_{8}=0.71\sigma$ when comparing the two HMCode2020 analysis. We note
that for the HMCode2020 analysis, the $S_{8}$ error from the
$E_{n}(\theta=$0$\stackrel{{\>{}^{\prime}}}{{\textstyle.}}$5$)$ data vector is
$\sim 15\%$ higher than the $S_{8}$ error from
$E_{n}(\theta=2.^{{}^{\prime}}0)$ data vector. This may be caused by chain-to-
chain variance noise. Defining the $S_{8}$ offset for the two HMCode2020
analyses in terms of the error on $E_{n}(\theta=2.^{{}^{\prime}}0)$, we find
$\Delta S_{8}=0.81\sigma$.. In this appendix we focus on understanding the
likely origin of a $0.7-0.8\sigma$ offset associated with scale cuts.

Figure 19: Estimating the expected bias $\Delta S_{8}$ when analysing a
Universe with no baryon feedback (upper right quadrant) or a Universe with
strong AGN feedback, modelled with Cosmo-OWLS:8.7 (lower left quadrant), with
a fiducial theoretical model that adopts HMCode2020 with $\log_{10}(T_{\rm
AGN}/{\rm K})=7.8$. We determine the signal-to-noise $\nu$ (Equation 136) for
three mock cosmic shear observables: DES Y3 $\xi_{\pm}(\theta)$ including
scale cuts (green), KiDS-1000 COSEBIs $E_{n}$ with $\theta_{\rm min}=$
2$\stackrel{{\>{}^{\prime}}}{{\textstyle.}}$0 (pink) and KiDS-1000 $E_{n}$
with $\theta_{\rm min}=$ 0$\stackrel{{\>{}^{\prime}}}{{\textstyle.}}$5 (blue).
This is calculated for the fiducial model over a range of $S_{8}$ values
(solid lines), and also for the trial cases of no feedback, and strong
feedback (data points). As expected the two scale cut statistics
($\xi_{\pm}(\theta)$ \- circle, $E_{n},\theta_{\rm min}=$
2$\stackrel{{\>{}^{\prime}}}{{\textstyle.}}$0 - star) produce the least biased
estimate when the underlying feedback differs from the model. The offset
between the two COSEBIs statistics is within the $1\sigma$ error for a KiDS-
like analysis (shaded).

In the interest of conducting a fast approximate analysis, we choose to
compress the tomographic cosmic shear data vector into a single signal-to-
noise estimate $\nu$ (Equation 136), following Appendix E.1. By noting that
the amplitude of the cosmic shear signal scales roughly as $\nu\propto
S_{8}^{2}$ (Jain & Seljak, 1997), we can write that the ratio $\nu_{\rm
A}/\nu_{\rm B}\propto\Delta S_{8}$, for the same noise realisation of
cosmology A and B whose $S_{8}$ values differ by a small quantity $\Delta
S_{8}=S_{8}^{\rm A}-S_{8}^{\rm B}$.

In Figure 19 we calibrate this relationship by calculating the ratio of
$\nu_{\rm test}/\nu_{\rm fiducial}$ for a series of mock KiDS-1000 and DES Y3
data vectors, using test cosmologies that differ from our fiducial
cosmological parameter set in Table 6 with $S_{8}^{\rm test}=S_{8}^{\rm
fiducial}+\Delta S_{8}$. We find a slightly different relationship for the DES
Y3-like $\xi_{\pm}(\theta)$ data vector (green) and the KiDS-1000-like $E_{n}$
data vectors (pink, blue), but a similar response for the two COSEBIs
statistic versions with $\theta_{\rm min}=$
0$\stackrel{{\>{}^{\prime}}}{{\textstyle.}}$5 (blue) and $\theta_{\rm min}=$
2$\stackrel{{\>{}^{\prime}}}{{\textstyle.}}$0 (pink).

To address the question, what offset in $S_{8}$ would we expect to measure
from the different cosmic shear statistics if our baryon feedback model was
incorrect, we calculate $\nu_{\rm test}/\nu_{\rm fiducial}$ for two further
tests. In both tests the cosmological model is fixed to the fiducial parameter
set and the baryon feedback model is altered. First we adopt a dark matter
only matter power spectrum with no baryon feedback suppression. The fiducial
mock simulates baryon feedback with $\log_{10}(T_{\rm AGN}/{\rm K})=7.8$,
reducing power relative to the dark matter only case by a factor of $\sim 5\%$
at $k=1h{\rm Mpc}^{-1}$. In this case $\nu_{\rm test}/\nu_{\rm fiducial}>1$
(see the upper right quadrant in Figure 19). For the second test we adopt a
model from the most extreme baryon feedback simulation to date Cosmo-OWLS:8.7
(Le Brun et al., 2014), which reduces power relative to the dark matter only
case by a factor of $\sim 20\%$ at $k=1h{\rm Mpc}^{-1}$, roughly a factor of
four times the level of suppression relative to the fiducial $\log_{10}(T_{\rm
AGN}/{\rm K})=7.8$ case. In this test $\nu_{\rm test}/\nu_{\rm fiducial}<1$
(see the lower left quadrant in Figure 19). Using the $\nu_{\rm A}/\nu_{\rm
B}\propto\Delta S_{8}$ relationship calibrated for each statistic in Figure
19, we can read off an estimate for $\Delta S_{8}$ for these two test cases.
As expected we find the COSEBIs $\theta_{\rm min}=$
0$\stackrel{{\>{}^{\prime}}}{{\textstyle.}}$5 statistic (blue) to be the most
sensitive to baryon feedback (the largest $\Delta S_{8}$ values), followed by
COSEBIs $\theta_{\rm min}=$ 2$\stackrel{{\>{}^{\prime}}}{{\textstyle.}}$0
(pink) and then $\xi_{\pm}(\theta)$ with DES Y3 scale cuts applied. The
difference seen between the COSEBIs $\theta_{\rm min}=$
2$\stackrel{{\>{}^{\prime}}}{{\textstyle.}}$0 and $\xi_{\pm}(\theta)$ scale
cut analyses demonstrates that the COSEBIs-defined scale cut from Appendix B
renders the cosmic shear analysis more sensitive to the $k$-scales impacted by
baryon feedback, compared to the Krause et al. (2021)
$\xi_{\pm}(\theta)$-defined scale cuts. It may be beneficial in the future to
develop the COSEBIs software to allow for redshift-dependent $\theta_{\rm
min}$-limits in order to further reduce the impact of baryon feedback without
incurring a significant loss of constraining power.

What is of interest for this appendix is the $S_{8}$ offset found between the
different statistics, rather than the offset relative to the fiducial input
value. If the underlying truth was a Universe with no baryon feedback, in a
noise-free analysis we would expect our best fit value for $S_{8}$ to decrease
by $\Delta S_{8}=0.012$ when changing from $\theta_{\rm min}=$
0$\stackrel{{\>{}^{\prime}}}{{\textstyle.}}$5 to $\theta_{\rm min}=$
2$\stackrel{{\>{}^{\prime}}}{{\textstyle.}}$0\. This corresponds to a $\sim
0.6\sigma$ decrease between the different KiDS-1000 constraints. In the data
we find a higher $S_{8}$ value when scale cuts are included, which is more
similar to the case where the underlying truth is Cosmo-OWLS:8.7. Here we
would expect to measure an increase of $\Delta S_{8}=0.008$ between the two
noise-free COSEBIs analyses, corresponding to a $0.4\sigma$ offset between the
different KiDS-1000 constraints.

It is important to recognise that these estimates are only approximate as they
compress the tomographic data vector into a single noise-averaged value,
erasing second-order information on scale and redshift dependence. In one
respect, however, they are conservative, as they do not account for the
marginalisation step over a wide baryon feedback prior in the KiDS-like
analysis that allows for models as strong as $\log_{10}(T_{\rm AGN}/{\rm
K})=8.0$. From this analysis we draw the conclusion that the $S_{8}$ offset we
measure between the KiDS-1000 COSEBIs $\theta_{\rm min}=$
0$\stackrel{{\>{}^{\prime}}}{{\textstyle.}}$5 and $\theta_{\rm min}=$
2$\stackrel{{\>{}^{\prime}}}{{\textstyle.}}$0 $E_{n}$ constraints could be
attributed to unmodeled extreme baryon feedback. This hypothesis only accounts
for roughly half of the observed difference, however, when considering Cosmo-
OWLS:8.7 as the ‘worst-case scenario’ for baryon feedback.

Asgari et al. (2021) report cosmological parameter constraints from a range of
two-point statistics measured from the same data set: $E_{n}(\theta_{\rm
min}=$0$\stackrel{{\>{}^{\prime}}}{{\textstyle.}}$5$)$, $\xi_{\pm}(\theta)$
and $C^{\rm EE}_{\epsilon\epsilon}(\ell)$. Using a series of 5000 noisy mock
simulations of KiDS-1000 created with SALMO727272SALMO:
https://github.com/Linc-tw/salmo. For information see section 4 of Joachimi et
al. (2021)., they quantify the cross-correlation between the different
statistics. The range of positive, zero and negative cross-correlations for
different components result from the changing sensitivity of each statistic to
a range of noisy $\ell$-scales. Asgari et al. (2021) find that they expect to
measure a $>1\sigma$ offset in $S_{8}$ constraints between the different
KiDS-1000 two-point statistics, 15% of the time (see similar findings for
analyses of HSC-Y1 mocks in Hamana et al., 2022). The $0.3\sigma$ offset
between the fiducial $\xi_{\pm}(\theta)$ and $E_{n}$ constraints seen in Table
13 is therefore expected.

We revisit the Asgari et al. (2021) analysis to quantify the expected offset
in $S_{8}$ between two COSEBIs analyses of the same data set using
$\theta_{\rm min}=$ 0$\stackrel{{\>{}^{\prime}}}{{\textstyle.}}$5 and
$\theta_{\rm min}=$ 2$\stackrel{{\>{}^{\prime}}}{{\textstyle.}}$0\. Using the
original set of SALMO mocks, we were able to determine the cross-correlation
between the two COSEBIs versions at essentially no CPU cost. We then took a
short-cut compared to the original analysis, using the $\nu_{\rm A}/\nu_{\rm
B}\propto\Delta S_{8}$ relationship derived in Figure 19 to determine the
expected offset in $S_{8}$ between pairs of $E_{n}$ statistics measured with
and without scale cuts in each SALMO mock. This differs from the CPU-intensive
approach of Asgari et al. (2021) where full CosmoSIS analyses were performed
for each mock and each statistic. By applying the same short-cut technique to
the set of statistics tested in Asgari et al. (2021), however, we have
verified that this approximate approach results in a similar distribution in
$\Delta S_{8}$ values. We find that a $\geq 1\sigma$ offset in $S_{8}$ between
the two COSEBIs measurements is expected 9% of the time. As such the two
COSEBIs versions are more correlated with each other than with
$\xi_{\pm}(\theta)$ and $C^{\rm EE}_{\epsilon\epsilon}(\ell)$. For the
$0.7\sigma$ difference that we observe in the HMCode2016 analysis, we find
that an offset of $\geq 0.7\sigma$ is expected 23% of the time. An offset of
$\geq 0.8\sigma$ is expected 14% of the time. As such the offsets we find
between the scale cut and original version of COSEBIs with KiDS-1000 can be
attributed to random noise fluctuations.

To the reader most familiar with the $\xi_{\pm}(\theta)$ statistic, this
result may seem counter-intuitive, as the scale cut that has been applied
removes only one data point from each tomographic bin. As we can see in Table
13, a $\xi_{\pm}(\theta)$ analysis with a $\theta_{\rm min}=$
2$\stackrel{{\>{}^{\prime}}}{{\textstyle.}}$0 scale cut applied indeeds makes
little impact, with $S_{8}$ changing by $\sim 0.2\sigma$. In contrast, the
COSEBIs statistic is calculated as a weighted average of finely binned noisy
measurements of $\xi_{\pm}(\theta)$ summed over the range $\theta_{\rm
min}<\theta<\theta_{\rm max}$. The shape of the oscillatory weight functions
depends on the angular range limits (see equations 28 and 38 of Schneider et
al., 2010). For $\theta_{\rm min}=$
0$\stackrel{{\>{}^{\prime}}}{{\textstyle.}}$5 and $\theta_{\rm min}=$
2$\stackrel{{\>{}^{\prime}}}{{\textstyle.}}$0 it is worth noting that the
COSEBIs weight functions differ significantly for $\theta\lesssim$
40$\stackrel{{\>{}^{\prime}}}{{\textstyle.}}$0.

##  Appendix G Constraints on all parameters, additional tables and figures

Summary: In this appendix we present additional results from the DES Y3 +
KiDS-1000 cosmic shear analysis to complement the primary results in Section
3.

Survey | Pipeline | $\chi^{2}_{\rm min}$ | $N_{\Theta}$ |  $p(\chi_{\rm min}^{2}|\nu_{\rm eff}$)
---|---|---|---|---
DES Y3 (Full area) | DES-like | 283.6 | 6.7 | 0.222
DES Y3 (Full area) | KiDS-like | 284.3 | 5.0 | 0.235
DES Y3 (Full area) | Hybrid | 284.2 | 5.4 | 0.231
DES Y3 (KiDS-excised) | DES-like | 285.3 | 8.0 | 0.187
DES Y3 (KiDS-excised) | KiDS-like | 286.5 | 5.1 | 0.207
DES Y3 (KiDS-excised) | Hybrid | 288.3 | 4.6 | 0.192
KiDS | DES-like | 83.7 | 8.3 | 0.078
KiDS | KiDS-like | 82.5 | 6.5 | 0.118
KiDS | Hybrid | 88.3 | 7.1 | 0.048
Joint | DES-like | 371.5 | 12.0 | 0.089
Joint | KiDS-like | 381.3 | 7.8 | 0.062
Joint | Hybrid | 378.0 | 9.6 | 0.068
Joint-1IA | DES-like | 377.2 | 10.4 | 0.068
Joint-1IA | KiDS-like | 373.9 | 8.8 | 0.094
Joint-1IA | Hybrid | 382.2 | 8.0 | 0.057
Table 14: Goodness of fit statistics listing the $\chi^{2}_{\rm min}$, the effective number of parameters, $N_{\Theta}$, and the goodness of fit statistic $p(\chi^{2}>\chi^{2}_{\rm min}|\nu)$ (see Section 2.7 for details). The Hybrid pipeline results can be compared to the DES-like and KiDS-like analysis of DES, KiDS and the joint survey. Parameter | DES Y3+KiDS-1000 | DES Y3 | KiDS-1000
---|---|---|---
$S_{8}$ | $0.790^{+0.018}_{-0.014}$ | $0.802^{+0.023}_{-0.019}$ | $0.763^{+0.031}_{-0.023}$
$\Omega_{\rm m}$ | $0.280^{+0.037}_{-0.046}$ | $0.297^{+0.040}_{-0.060}$ | $0.270^{+0.056}_{-0.102}*$
$\sigma_{8}$ | $0.825^{+0.067}_{-0.073}$ | $0.816^{+0.076}_{-0.085}$ | $0.833^{+0.133}_{-0.146}$
$h$ | $0.746^{+0.069}_{-0.029}*$ | $0.738^{+0.072}_{-0.040}*$ | $0.749^{+0.068}_{-0.029}*$
$n_{\rm s}$ | $0.973^{+0.074}_{-0.072}*$ | $0.979^{+0.093}_{-0.062}*$ | $0.984^{+0.102}_{-0.053}*$
$\Sigma m_{\nu}$ | $0.326^{+0.110}_{-0.187}*$ | $0.317^{+0.139}_{-0.223}*$ | $0.327^{+0.152}_{-0.215}*$
$\omega_{\rm b}$ | $0.022^{+0.002}_{-0.003}*$ | $0.022^{+0.002}_{-0.003}*$ | $0.022^{+0.002}_{-0.003}*$
$\mathrm{log}_{10}(T_{\mathrm{AGN}}/{\rm K)}$ | $7.627^{+0.152}_{-0.299}*$ | $7.632^{+0.170}_{-0.289}*$ | $7.623^{+0.153}_{-0.293}*$
$A_{\rm IA}^{\mathrm{DES}}$ | $-0.025^{+0.576}_{-0.285}$ | $0.323^{+0.432}_{-0.369}$ | -
$\eta^{\mathrm{DES}}$ | $2.047^{+2.953}_{-0.832}*$ | $1.774^{+3.166}_{-0.940}*$ | -
$A_{\rm IA}^{\mathrm{KiDS}}$ | $1.038^{+0.543}_{-0.516}$ | - | $0.650^{+0.884}_{-0.539}$
$\eta^{\mathrm{KiDS}}$ | $2.211^{+2.588}_{-0.958}*$ | - | $2.022^{+2.937}_{-0.868}*$
$\Delta z_{1}^{\mathrm{DES}}$ | $0.001^{+0.015}_{-0.015}*$ | $-0.004^{+0.016}_{-0.015}*$ | -
$\Delta z_{2}^{\mathrm{DES}}$ | $0.014^{+0.012}_{-0.011}*$ | $0.011^{+0.011}_{-0.011}*$ | -
$\Delta z_{3}^{\mathrm{DES}}$ | $-0.008^{+0.010}_{-0.010}*$ | $-0.004^{+0.010}_{-0.010}*$ | -
$\Delta z_{4}^{\mathrm{DES}}$ | $-0.006^{+0.015}_{-0.014}*$ | $0.000^{+0.014}_{-0.014}*$ | -
$\Delta z_{1}^{\mathrm{KiDS}}$ | $-0.001^{+0.011}_{-0.010}*$ | - | $0.000^{+0.010}_{-0.010}*$
$\Delta z_{2}^{\mathrm{KiDS}}$ | $-0.010^{+0.010}_{-0.011}*$ | - | $-0.008^{+0.011}_{-0.012}*$
$\Delta z_{3}^{\mathrm{KiDS}}$ | $0.016^{+0.010}_{-0.011}*$ | - | $0.016^{+0.011}_{-0.010}*$
$\Delta z_{4}^{\mathrm{KiDS}}$ | $0.017^{+0.008}_{-0.009}*$ | - | $0.015^{+0.008}_{-0.008}*$
$\Delta z_{5}^{\mathrm{KiDS}}$ | $-0.006^{+0.009}_{-0.010}*$ | - | $-0.007^{+0.009}_{-0.010}*$
$m_{1}^{\mathrm{DES}}$ | $-0.007^{+0.009}_{-0.009}*$ | $-0.007^{+0.009}_{-0.009}*$ | -
$m_{2}^{\mathrm{DES}}$ | $-0.021^{+0.008}_{-0.007}*$ | $-0.020^{+0.008}_{-0.008}*$ | -
$m_{3}^{\mathrm{DES}}$ | $-0.022^{+0.007}_{-0.007}*$ | $-0.024^{+0.007}_{-0.008}*$ | -
$m_{4}^{\mathrm{DES}}$ | $-0.035^{+0.007}_{-0.007}*$ | $-0.037^{+0.007}_{-0.007}*$ | -
Table 15: Mean marginal constraints with 68% credible intervals for the
$\Lambda$CDM parameters (upper), astrophysical nuisance parameters (middle)
and data calibration parameters (lower). Cosmic shear results are reported
from the joint DES Y3 and KiDS-1000 (left), DES Y3 (middle) and KiDS-1000
(right) analyses. Parameters marked ∗, are considered to be unconstrained by
the data as their marginal posterior, at either prior edge, exceeds 13% of the
peak posterior probability.

  * •

Table 14 presents the goodness of fit statistics for the DES-like and KiDS-
like analyses, listing the $\chi^{2}_{\rm min}$, the estimated effective
number of parameters, $N_{\Theta}$, and the goodness of fit statistic
$p(\chi^{2}>\chi^{2}_{\rm min}|\nu)$ for a range of different analyses of the
DES, KiDS and joint-survey data vectors, with $N_{\rm data}=273,75,348$
respectively. We remind the reader that the Joachimi et al. (2021) approach
used to estimate $N_{\Theta}$ is subject to noise at the level of $\sim
0.2N_{\Theta}$.

  * •

In Figure 20 we present the joint-survey constraints across the full
cosmological parameter space, comparing the fiducial Hybrid NLA-z analysis
with the alternative TATT analysis. We also include the marginal posterior for
the baryon feedback parameter $T_{\rm AGN}$. Only $S_{8}$, $\Omega_{\rm m}$
and $\sigma_{8}$ are constrained by the data. The axis limits on $h,n_{\rm
s},\Sigma m_{\nu}$, $\omega_{\rm b}$ and $T_{\rm AGN}$ are set by the priors.
As we have chosen to sample in $S_{8}$, the prior space is straightforward to
visualise as rectangles in each 2D marginal bounded by the prior edges.
Throughout the paper we use GetDist737373GetDist:
https://getdist.readthedocs.io (Lewis, 2019) to plot quantities of interest
from our chains of Monte Carlo samples. GetDist applies a linear boundary
kernel and multiplicative bias correction747474We use the GetDist settings:
boundary_correction_order:0, mult_bias_correction_order:1. to remove the bias
introduced when the chain-smoothing kernel passes over a prior boundary. When
the hard boundaries from the priors are aligned with the parameter
coordinates, as in our case, the boundaries can be defined using the
GetDist.MCSamples ‘ranges’ option. With the prior ranges set, the GetDist
contoured credible intervals do not artificially close for prior-informed
parameters.

  * •

In Figure 21 we present the joint-survey constraints across the astrophysical
nuisance parameter space for the alternative TATT analysis. Only $S_{8}$,
$a_{1}$ and $a_{2}$ are constrained by the data. We can see that the
preference for a lower $S_{8}$ value in the TATT analysis, relative to our
fiducial analysis, is primarily driven by the IA model allowing for a
significant tidal torque amplitude, the $a_{2}$ parameter, when there is
strong redshift evolution with high values for $\eta_{2}$. We find a similar
preference for strong redshift evolution of the IA model in the fiducial NLA-z
analysis (see Figure 3), but in contrast to the TATT analysis, the $\eta_{\rm
IA}$ amplitude is largely uncorrelated with $S_{8}$. We see similar behaviour
for the IA model parameters for the DES (pink) and KiDS (black) surveys. As
with the fiducial analysis, KiDS prefers higher values for $A_{1}$ compared to
DES.

  * •

Table 15 reports the mean 1D marginal values for each parameter for the
fiducial Hybrid analyses of DES Y3 and KiDS-1000. We use ∗ to mark parameters
that are unconstrained by the data. In these cases the tabulated constraint
comes from the prior, as the marginal posterior probability, at either prior
edge, exceeds 13% of the peak posterior value (see appendix A of Asgari et
al., 2021). We constrain $S_{8}$, $\sigma_{8}$, and the intrinsic alignment
amplitude $A_{\rm IA}$ for each survey, and $\Omega_{\rm m}$ for the joint
survey and the DES-only analysis. We have not presented figures of the
marginal posteriors for the shear and redshift calibration nuisance parameters
$m$ and $\Delta z$, as we recover the peak and width of the Gaussian prior for
the majority of these parameters (compare the values757575In this paper we
define $\Delta z=\langle n^{\rm estimate}(z)\rangle-\langle n^{\rm
true}(z)\rangle$ which we note has a different sign convention to the CosmoSIS
parameter DELTA_Z_OUT$=-\Delta z$. in Table 1 and Table 15). A notable
exception is the redshift calibration parameter for the second tomographic
bin, $\Delta z_{2}$, for both the DES and KiDS data. In the DES-only, KiDS-
only and joint-survey Hybrid analysis the mean $\Delta z_{2}$ of each survey
is $\pm 1\sigma$ offset from the Gaussian prior peak. This is consistent with
the previous findings of both survey teams (see for example figure 17 of Amon
et al., 2022). Both teams show consistent $S_{8}$ results when excluding
individual tomographic bins from the fiducial analysis (see also figures 16
and 18 in Amon et al., 2022, for the stability of the intrinsic alignment
parameter constraints when removing redshift bins and varying the redshift
calibration methodology). Using the Köhlinger et al. (2019) ‘three tier’ set
of internal consistency statistics, Asgari et al. (2021) flag the second
KiDS-1000 tomographic bin as an outlier (see their appendix B.2). As this bin
carries insignificant levels of cosmological information, however, the results
are unchanged by the inclusion or exclusion of this tomographic bin.

Figure 20: Joint survey posteriors for the $\Lambda$CDM cosmological
parameters and the HMCode2020 $T_{\rm AGN}$ parameter for baryonic feedback.
We compare the fiducial Hybrid (NLA-z) analysis in pink to a TATT analysis
with black dashed lines. The inner 68% and outer 95% credible intervals are
shown contoured. In the case of the parameters $h,n_{\rm s},\Sigma m_{\nu}$,
$\omega_{\rm b}$ and $T_{\rm AGN}$, the axis limits are set by the prior
boundaries. Figure 21: Joint survey posteriors, comparing constraints on the
TATT intrinsic alignment model (see Section 2.2) from the DES Y3 (pink) and
KiDS-1000 (black-dashed) surveys with the inner 68% and outer 95% credible
intervals contoured. In all but the $S_{8}$ parameter, the axis limits are set
by the prior boundaries.