# Understanding the Kerala Floods of 2018:
Role of Mixed Rossby-Gravity Waves

Kiran S. R.
Lecturer in Civil Engineering, Central Polytechnic College
Department of Technical Education
Thiruvananthapuram - 695013
Kerala, India

_Submitted on $17^{th}$ March 2021_
Corresponding author<EMAIL_ADDRESS>ce_kiran@cpt.ac.in)_

###### Abstract

Kerala, the south-west coastal state of India, was ravaged by a series of
floods during the South-West Monsoon of 2018. The season was marked by
severely anomalous rainfall trends, with upto 150 mm of departures from the
mean daily precipitation in the northern districts of the State. Although,
there were many studies about the hydrogeological factors which aggravated the
floods in Kerala, no attempt was made to delve into the physics which actually
resulted in anomalous precipitation during the year. This study intends to
document the dynamical phenomenon which caused the Kerala Floods of 2018. The
westward propagating convectively-coupled Mixed Rossby-Gravity (MRG) waves
were excited by the synoptic disturbances of the tropical Pacific at the
pressure level of 700 hPa, during the Indian Monsoon of 2018. They travelled
across the Indian Ocean in two significant modes – a predominant slow moving
wave of 20-40 days period (as Madden-Julian Oscillation) and a secondary
faster wave of 5-8 days period. They are characterised by vertical phase
propagation to the upper troposphere, a precursor to deep convection and
intense precipitation. Further, the propagation of these waves through a
medium enhances its relative vorticity and the gyres or circulations thus
formed are symmetrical about the equator. Consequently, the meanders in the
wind field and widening of the Intertropical Convergence Zone were observed.
The MRG waves, especially the slow mode induced divergence in the wind field,
which fueled convection in tropics and brought very heavy rainfall to the
State of Kerala in 2018.

_Keywords_ Kerala $\cdot$ Floods $\cdot$ Monsoon $\cdot$ Intertropical
Convergence Zone $\cdot$ Mixed Rossby-Gravity Waves

##  1 Introduction

Climate of any place is the product of inextricable intertwining of the
atmosphere, ocean, land, ice and biology. Climate change has become a buzzword
today and is marked by frequent occurrences of extreme weather events – severe
floods and droughts (Jentsch and Beierkuhnlein, 2008; Bergquist et al., 2019).
It is observed that the frequency of occurrences of short-duration extreme
precipitation events have increased in different parts of the world, including
India.

Kerala, the State on south-west coast of India, experienced torrential
rainfall and devastating floods in the year 2018 during Boreal Summer. The
South-West Monsoon season observed an exceptionally high precipitation in the
State with over 15%, 18% and 164% above normal rainfall during the months of
June, July and August respectively (IMD, 2018). Kerala received a total of 821
mm of precipitation only in the month of August, which is the highest recorded
rainfall for the month after 1931. A series of unprecedented deluge events
ensued the heavy downpour and 13 out of the 14 districts in the State
witnessed floods of disastrous proportions. Further, the State witnessed
recurrent flood events, during the same season, in the subsequent years 2019
and 2020.

Mishra and Shah (2019) observed that there is no increase in mean and extreme
precipitation in Kerala over the last six decades and the extreme rainfall
event during August 2018 is driven by anomalous atmospheric conditions due to
climate variability rather than anthropogenic factors. They concluded that the
severity of the Kerala flood of 2018 was exacerbated by several other factors
such as land use/ land cover change, antecedent hydrologic conditions,
reservoir storage and operations, encroachment of flood plains, etc. Lal et
al. (2020) attempted to perform inundation mapping of Kerala, to identify the
flood-affected regions in the State. Sudheer et al. (2019) analyzed the role
of dams and reservoir operations, on the flood of August 2018. They deduced
that the role of water, released from the major reservoirs in the Periyar
River Basin, to have resulted in the flood havoc was minimal. Hunt and Menon
(2020), using Weather Research and Forecasting (WRF) model in tandem with a
hydrological model, performed a broad analysis to underscore climate change
and concluded that the precipitation may be 18% heavier for the pre-industrial
period (due to recent weakening of monsoon low-pressure systems), and 36%
heavier for a futuristic climate of 2100 (due to moistening of the tropical
troposphere).

Previous studies, about the Kerala Floods of 2018, have neither documented the
observations nor attempted to analyze deeper into the physical phenomenon
responsible for the disaster. Hence, any investigation into the anomalous
climate dynamics, which triggered the deadly flood events in the State, is yet
to pull off. Further, the recurrence of severe flood events in Kerala in the
recent years arose the need to identify this anomaly at the earliest, which
became the prime motivation for this study.

##  2 Data & Methodology

The daily rainfall data over the Indian Ocean region is obtained from Tropical
Rainfall Measuring Mission (TRMM) at a spatial resolution of $0.25^{o}$
$\times$ $0.25^{o}$, for the period from 2011 to 2018. It gives the daily
accumulated precipitation (combined Microwave-IR) estimate with gauge
calibration performed over land in millimetres. The globally gridded data of
zonal and meridional winds at a spatial resolution of $2.5^{o}$ $\times$
$2.5^{o}$ from 2011 to 2018 is obtained from NCEP/NCAR Reanalysis 1, which
uses data assimilation to incorporate observations and numerical weather
prediction (NWP) model output, for accurate representation of the state of
Earth’s atmosphere. This gives horizontal wind (zonal and meridional) data (in
metres per second) for different vertical pressure levels such as 1000, 925,
850, 700, 600, 500, 400, 300, 250 and 200 in hectopascals. In order to analyse
the probable effects of convection, omega or temporal derivative of pressure
(in pascals per second), which is an estimate of vertical velocity for an
atmosphere in pressure coordinate system, is also obtained for different
vertical levels, as mentioned above, from NCEP/NCAR Reanalysis.

In this study, the daily departures from the mean of atmospheric variables
during 2018 are determined to quantify the gravity of anomalous precipitation
in Kerala. This is estimated by deducting the average of daily data of the
atmospheric variables from 2011 to 2017 (7year daily mean) from the daily data
for the year 2018. In order to compute the width of Intertropical Convergence
Zone (ITCZ), a threshold rainfall of 50mm/day is chosen and the latitudinal
extend over which the spatial continuous rainfall outstrip the threshold along
the desired meridian is deemed to be the apparent width of ITCZ. Computation
of the relative vorticity involves determination of spatial derivatives of
winds, which is performed by means of simple finite-backward differencing
method, as executed by Kiran (2017) to determine the curl of wind stress in
Andaman Sea. Hovmöller diagrams, which depict the temporal variation of
atmospheric parameters in the zonal or meridional direction, are used to
identify and explain atmospheric waves and its characteristics. Fast Fourier
Transform (FFT) is performed on time series of atmospheric variables or its
anomalies for the ease of representation in frequency domain. The dominant
harmonics, thus identified, are filtered out using Band-Pass Filters (BPF) for
further analysis.

##  3 Results & Discussions

###  3.1 Rainfall over Kerala in 2018

The Kerala Flood of 2018 was the result of intense precipitation during the
South-West Monsoon over the Indian subcontinent during the months of June to
August. The daily rainfall during the period at Mankadavu (Idukki, Central
District of Kerala) and Pavannur (Kannur, Northern District of Kerala), as
obtained from TRMM dataset, is shown in Figure 1 (top). Pavannur leads in
terms of rainfall intensity compared to Mankadavu, but both the locations
captured the major trends in precipitation. It is observed that daily rainfall
in excess of 50mm, occurred at both the locations during the periods $8^{th}$
\- $11^{th}$ June, $3^{rd}$ \- $7^{th}$ July, $16^{th}$ \- $26^{th}$ July and
$4^{th}$ \- $12^{th}$ August, although Mankadavu missed the peak rainfall
during the first week of July. The above time durations are considered as
occurrences of extreme weather events and are designated as T1, T2, T3 and T4
respectively (refer Table 1 to compare the rainfall trends at both the
stations). A very heavy rainfall of 200 mm/day and above occurred by the end
of July and beginning of August, which aggravated the flood events in the
State. To quantify the rainfall anomalies, the precipitation data from TRMM
for 7 years (2011 to 2017) is daily-averaged and subtracted from the rainfall
of 2018 (Fig. 1 (bottom) & Table 1). The year 2018 saw an extremely anomalous
monsoon rainfall during T1, T2, T3 and T4, with excess daily precipitation as
high as 150 mm, and hence, this phenomenon requires to be investigated.

Figure 2 shows the spatial distribution of TRMM precipitation on the days of
peak rainfall each during T1, T2, T3 and T4 respectively, overlaid with
NCEP/NCAR Reanalysis surface wind vectors. On all these days, the entire
extend of the state of Kerala received abnormally high precipitation
(exceeding 50 mm per day) compared to its neighborhood. It may be observed
that the rainfall intensified on the windward side of Western Ghats, with the
wind vectors oriented mostly normal to the chain of mountains. This
underscores the role of orography in galvanizing Monsoon precipitation. On 5th
July 2018, the heavy rainfall is spread only across northern part of the
State, which explains why the anomaly during T2 was negligible at Mankadavu.
Further, one observes very severe horizontal wind shear and the associated
tendency for cyclonic circulation over the region on the days of peak
rainfall, especially during T3 and T4. Here, the positive relative vorticity
sets high prospects for good rainfall.

###  3.2 Analysing the behaviour of ITCZ

Intertropical Convergence Zone (ITCZ) or the Equatorial trough is an
atmospheric system responsible for the large-scale rainfall over tropical
oceans. It functions as the ascending limb of the meridional circulation
(known as Hadley cell) at lower latitudes of the atmosphere. Indian monsoon is
the manifestation of seasonal migration of ITCZ. Further, it is associated
with a series of northward propagation of the cloud bands from the equatorial
Indian Ocean onto the mainland (Riehl, 1954; Charney, 1969; Gadgil, 2018).

The location of ITCZ for the given year is chosen as the latitude
corresponding to the meridional maximum daily precipitation (Byrne et al.,
2018). The daily variation of the position of ITCZ along $76^{o}$E, the
meridian passing through the State of Kerala, is shown in Figure 3 (a), in
comparison with the 7 year (2011-2017) daily mean location of the Equatorial
trough. During T1, T2, T3 and T4, the ITCZ is positioned over Kerala, although
it considerably varies in location and even fell out-of-phase with the mean
during T3 and T4. Figure 3 (b) is the evolution of meridional rainfall maxima,
i.e., the peak rainfall within ITCZ in 2018. One finds the ITCZ rainfall
itself is the greatest during T1, T2, T3 and T4, as against the mean. Further,
there was a simultaneous increase in the width of ITCZ during T3 and T4 (Fig.
3 (c)), where the trough width is obtained as the latitudinal extend about the
precipitation maxima along which the daily rainfall exceeds a threshold of
50mm. Here, the widening of ITCZ may be attributed to the perturbations caused
to the south-westerly winds, by any local or remote phenomenon (Byrne and
Thomasi, 2019). Hence, ITCZ was approximately positioned over Kerala during
the occurrences of extreme rainfall, which widened and delivered heavy
precipitation over the State.

###  3.3 Understanding the waves and its dynamics

Indian Monsoon is said to be strong if the zonal component of the Somali Jet
is greater by more than a standard deviation of the long-term mean. The zonal
component of the wind is generally stronger during the season, analysing which
shall render a better insight into the associated intraseasonal variability
(Gadgil, 2018; Murtugudde et al., 2007). NCEP/NCAR zonal wind anomalies over
Kerala (averaged over $75^{o}$E to $77^{o}$E & $8^{o}$N to $12^{o}$N) in 2018,
is shown in Figure 4. In comparison to the surface, the westerly winds
intensify during T1, T2, T3 and T4 at 800 to 600 hPa pressure levels, while
the negative anomalies dominate 400 hPa to 700 hPa during the early weeks of
June and July. The occurrence of high positive zonal wind anomalies as high as
8 m/s at an altitude of average pressure 700hPa is concomitant with the
periods of anomalous precipitation, i.e., during T1, T2, T3 and T4. The
perturbations tend to propagate vertically upwards by the end of July from 700
hPa to upper troposphere and then downwards to the surface. This is
characteristic of the vertical propagation of waves in atmosphere. Since, the
wind anomalies appear to concentrate near the pressure level of 700 hPa, it
shall be chosen as a reference for further investigation.

####  3.3.1 Source of perturbation

Figure 5 gives a neat picture of the evolution of NCEP/NCAR zonal wind
anomalies of 2018 at 700 hPa along $10^{o}$N latitude, which passes through
Kerala. This Hovemöller diagram exposes westward propagating synoptic
disturbances or equatorial waves, which originated from East Atlantic and
Pacific Ocean basins, during the months of June to September. They have
wavelengths of the range 2500 - 5000 km, and manifest as shallow disturbances
(as large wind anomalies at 700-850 hPa levels) embedded in easterly trade
winds. These may form as a result of instabilities in the mean flow or due to
barotropic instability of the ITCZ (Reed and Recker, 1971; Ferreira and
Schubert, 1997; Serra et al., 2008), yet the dynamics of initiation of
tropical synoptic disturbances hitherto remains a mystery. These westward-
propagating signals travelled to eastern end of Equatorial Indian Ocean and
perturbed the wind system there. As a result, the zonal wind anomalies over
the Indian Ocean basin ($50^{o}$E to $100^{o}$E) intensified.

####  3.3.2 Mixed Rossby-Gravity (MRG) Waves

The enlarged segment of the Hovemöller diagram (Fig. 5) exposes the
intensification of zonal wind anomalies accompanied by westward propagation of
signals between 60E and 120E. This reveals synoptic scale waves which
originate near $150^{o}$E (western Equatorial Pacific) and propagate westwards
till $60^{o}$E (near the coast of Oman), via Bay of Bengal and Arabian Sea.
These propagations seem to amplify and coincide exactly with the occurrences
of anomalous rainfall in Kerala, especially during T3 and T4. From Figure 6,
the NCEP/NCAR meridional wind anomalies over the tropics ($50^{o}$E to
$150^{o}$E) support the same events of waves (phase) propagating westward, and
the anomalies change sign eastwards (wave group). The phase propagates
westward with a speed ($C_{p}$) of 10 m/s, while the wave packets seem to
propagate eastward with a group velocity ($C_{g}$) of 6 m/s. The wave
propagation characteristics resemble Mixed Rossby-Gravity (MRG) waves, which
are westward-propagating equatorially-trapped waves with synoptic time scales.
Earlier studies identified lower-tropospheric MRG waves with 3–8 days period
in the central Pacific having wavelength of 6000–10000 km, a westward phase
speed of 10–20 m/s, and an eastward group speed of nearly 5 m/s (Dickinson and
Molinari, 2002; Raupp and Dias, 2005; Chen and Huang, 2009; Yang and Ingersol,
2011). Therefore, MRG waves, triggered by the westward propagating synoptic
disturbances of the Pacific, intensified the monsoon rainfall over Kerala as
they propagated across the equatorial Indian Ocean region in 2018.

####  3.3.3 Characteristics of MRG Waves

The NCEP/NCAR omega anomalies of 2018 averaged over $75^{o}$E to $77^{o}$E &
$8^{o}$N to $12^{o}$N (Fig. 7) show clear signs of deep convective activity
near the coast of Kerala in the Summer monsoon of 2018, especially during T3
and T4. It elicits an interesting wave character; there occurs upward phase
propagation (of negative omega or rising motion of air) from 700 hPa pressure
level to upper troposphere which supports deep convection and the consequent
intense precipitation is concurrent with the periods of occurrence of
anomalous rainfall. Besides, the wave group seem to propagate upward and
downward across the height of troposphere, which substantiates the
observations discussed on the anomalies of zonal velocity (Fig. 4). The rising
motion of the wave group appears to be a precursive event to the anomalous
precipitation, while the downward propagation follows it.

It is important to quantify the relative vorticity of the wind field, as waves
propagate through the medium. The strong shear in winds during the monsoon of
2018 was already discussed with Figure 2, which shall favour an increase in
the flow vorticity. The daily spatial variation of anomalies in the relative
vorticity at 700 hPa level overlaid with winds at the same pressure level from
$16^{th}$ July to $20^{th}$ July 2018 are shown in Figure 8. The winds near
Western Pacific possess very high vorticity of atleast 3 $\times$ $10^{-5}$
$s^{-1}$ on $16^{th}$ July, which tends to shift westward in the subsequent
days and intensify over the Bay of Bengal and Arabian Sea on $19^{th}$ July.
The wind field appears to possess a meandering or wavy pattern, which may be
attributed to the westward propagating MRG waves.

###  3.4 Spectral Analysis

The presence of prominent periodicities in the wind field may be determined by
Fourier Transforms, which converts time domain into frequency domain, as
discussed under Section 2. Fast Fourier Transform (FFT) of the relative
vorticity anomalies of 2018 for the region $60^{o}$E to $100^{o}$E & $5^{o}$N
to $20^{o}$N is performed and the spatially averaged resultant amplitude
spectrum is given in Figure 9. It may be observed that the dominant periods of
variability associated with relative vorticity belong to the ranges such as
5-8 days, 12-16 days and 20-40 days, in addition to the semiannual variability
(180 days period corresponds to the much acknowledged strong and sustained
westerly winds over the Equatorial Indian Ocean region during Spring and Fall
(Wyrtki, 1973)). Hence, the winds at 700 hPa pressure level are bandpass
filtered for the above ranges of periods, for further analysis (Fig. 10). Out
of the above harmonics, 12-16 days variability is seldom discussed here, owing
to its weak contribution to MRG wave dynamics. Therefore, the analysis of the
horizontal wind field at the equatorial Indian Ocean and West Pacific (over
$60^{o}$E to $120^{o}$E belt, averaged about latitudinal band $5^{o}$N to
$15^{o}$N) at 700 hPa, filtered separately for the bands 20-40 days and 5-8
days, are discussed here.

####  3.4.1 20-40 day mode

From the Hovemöller diagram (Fig. 10 (first & second figures)), one observes
westward propagating phase between $60^{o}$E and $90^{o}$E in both zonal and
meridional winds at 700 hPa level, with its amplitude maximized during T3 and
T4. The waves originate at the eastern end of the equatorial Indian Ocean as a
result of intensification of winds over the West Pacific region ($100^{o}$E to
$120^{o}$E) perturbed by synoptic scale disturbances, as discussed under
Section 3.3.2. Hence, the Mixed Rossby-Gravity (MRG) waves excited in the
Indian Ocean exist as a dominant mode of 20-40 days periodicity. These MRG
waves possess a westward phase speed of 3 m/s in zonal winds and 6 m/s in
meridional winds. These disturbances tend to propagate slowly in eastward
direction at about 1 m/s. It is interesting to note that the Madden–Julian
oscillation (MJO), a prominent convective disturbances in the tropics, exist
with the same time period. Moreover, Yang and Ingersol (2011) explained that
MJO is a wave packet formed as a narrow frequency band of MRG waves, and MJO
propagates eastward with the group velocity of MRG waves. Hence, MJO
manifested as slow-moving MRG waves of 20-40 day variability in the equatorial
Indian Ocean, especially during T3 and T4.

####  3.4.2 5-8 day mode

From the Hovemöller diagram (Fig. 10 (third & fourth figures)), MRG waves of
shorter period (5-8 days) is also observed to exist over the Equatorial Indian
Ocean. This secondary mode of MRG waves are faster (phase speed of around 12
m/s) but of relatively lower power density, compared to the primary mode of
MRG waves (20-40 days). The filtered zonal winds expose the secondary wave
mode during T3, while the meridional winds reveal the same during T4. Hence,
MRG waves existed in two significant modes in the tropical atmosphere during
Indian Monsoon of 2018 - a slow but dominant wave mode of 20-40 days period
and a fast mode of 5-8 days period.

The above discussion emphasize the prominence of Mixed Rossby-Gravity (MRG)
waves, which manifested as dominant modes of variability over Equatorial
Indian Ocean in 2018, during the periods of anomalous rainfall over the State
of Kerala. These equatorially-trapped waves are convectively-coupled and
possess a characteristic wave structure; these waves manifest as synoptic
scale gyres centered over as well as symmetric about the equator, but with
anti-symmetric convective tendencies (Kiladis et al., 2009). Spatial
distribution of relative vorticity filtered for periods 20-40 days and 5-8
days each (Fig. 11) show symmetric gyres or vortices positioned over equator,
on the days corresponding to peak rainfall during T1, T2, T3 and T4. It may be
observed that the circulation remained symmetrical between $15^{o}$N and
$15^{o}$S. For example, on $18^{th}$ July 2018 (during T3), a gyre of
clockwise circulation (shown in blue shade; it has anticyclonic vorticity in
the Northern Hemisphere and cyclonic vorticity in the Southern Hemisphere and
represents a wave crest) occupied the equator symmetrically at $85^{o}$E
accompanied by a cyclonic vortex (wave trough) over the Indian peninsula and
an anti-cyclonic vortex (wave trough) over Southern Indian Ocean. Similarly
oriented set of crest and troughs of the MRG wave may be observed on $8^{th}$
August 2018 (during T4) at $75^{o}$E. The characteristic length scale of these
vortex structures differ for each modes of MRG waves. The primary mode (20-40
days) of MRG wave has a significant influence over the wind field of 700 hPa
level compared to the secondary mode (5-8 days). As explained by Kiladis et
al. (2009), MRG wave modes are convectively-coupled and dilate the wind field
to favour rising and sinking of air parcels. This shall be investigated with
the help of Figure 12, which gives the spatial distribution of bandpass-
filtered omega anomalies at 700 hPa level for each MRG wave modes, and shall
now aid to assess the convective tendencies forced by MRG wave modes during
the days of peak anomalous rainfall. Figure 12 (top) reveals that the primary
mode brought the entire State of Kerala under a deep convective zone (high
negative omega anomalies) during T1, T3 and T4, while the secondary mode
caused convection only to the northern part of the State during T2. Upon
comparison of spatial distributed omega anomalies (Fig. 12) with that of the
precipitation (Fig. 2), it may be concluded that for the primary mode of MRG
waves, the omega anomalies bears high spatial correlation with the rainfall on
the days of peak anomalies in precipitation. Hence, the MRG wave mode of 20-40
days period had a very profound impact on the rainfall across Kerala in 2018.
Yet, the faster MRG wave (secondary mode) cannot be neglected, as it assumed
significance during T2, when only the northern part of the State fell within
the convective zone. Therefore, both the primary and secondary modes of MRG
waves in superposition, played a pivotal role in driving convection during the
period of anomalous rainfall in the State of Kerala in 2018.

##  4 Summary & Conclusions

The State of Kerala experienced a spate of flood events during the South-West
Monsoon of 2018. It is much more than a hydrogeological problem, especially in
the pretext of growing concerns of changing climate and its serious
ramifications. Hence, this study specifically unearths the dynamical
phenomenon responsible for the anomalous course of events.

Kerala received very high precipitation in 2018, especially during the fag end
of July (T3) and beginning of August (T4). The northern part of the state
received more anomalous rainfall, with daily precipitation departures as high
as 150 mm on $8^{th}$ August. Western Ghats played a major role to limit the
spatial extend of extreme precipitation to its windward side. The
Intertropical Convergence Zone (ITCZ) cloud band maxima shifted a large number
of times to the mainland and centered over the State, but underwent
substantial fluctuations in its position during T3 and T4. Besides, the ITCZ
widened in tune with high meridional rainfall maxima during periods of
anomalous precipitation. This interesting, yet unprecedented behaviour of
Equatorial trough motivates one to foray into the atmospheric dynamics which
perturbed the tropical convective system.

The synoptic-scale disturbances of the Pacific or Atlantic Ocean travelled
westward and perturbed the tropical dynamics to cause heavy downpour over
Kerala in 2018. They triggered the equatorially-trapped Mixed Rossby-Gravity
(MRG) Waves over the eastern Equatorial Indian Ocean at 700 hPa pressure
level, and propagated westward with an average phase speed of 10 m/s to the
strike coast of Africa. They possessed an eastward group velocity of around 6
m/s. These waves tend to escalate the relative vorticity as they propagate
through the medium and meander the wind field. The MRG waves perturbed the
ITCZ too; the westerly wave group caused the convective trough to widen.
Further, their phase propagates vertically upward to the upper troposphere
accompanied by intense updraught of air, thereby fueling convection in the
tropics.

The Equatorial Indian Ocean basin resonated to the two convectively-coupled
modes of MRG waves – a primary slow-moving wave of period 20-40 days and a
secondary faster wave of period 5-8 days. The slow westward propagating wave
(phase speed of 3-6 m/s) has the period identical to that of Madden-Julian
Oscillations (MJO), which is a characteristic MRG wave group travelling
eastward. The westward propagating faster waves are characteristic MRG waves
with phase speed of around 12-13 m/s. Both the wave modes set up vortices
(both cyclonic and anti-cyclonic) in the medium and are hence responsible for
meandering of winds. These circulations or atmospheric gyres position
symmetrically over the equator, as they propagate across the basin and dilate
the wind field to establish zones of convection in the tropics. The primary
mode of MRG waves (20-40 days period) is mainly responsible for stimulated
convection which resulted into large anomalies in rainfall over Kerala in
2018. The faster secondary mode, although of relatively low power density,
could significantly describe the spatial distribution of convective zones
during first week of July. Therefore, tropical atmosphere resonated to the two
MRG wave modes in superposition and further, catalyzed convection over the
equatorial Indian Ocean, which resulted in heavy precipitation over the State
of Kerala in 2018.

## References

  * Jentsch and Beierkuhnlein [2008] Anke Jentsch and Carl Beierkuhnlein.  Research frontiers in climate change: Effects of extreme meteorological events on ecosystems.  _Competes Rendus Geoscience, Vol. 340, Issues 9-10_ , 2008.
  * Bergquist et al. [2019] Magnus Bergquist, Andreas Nilsson, and P. Wesley Schultz.  Experiencing a severe weather event increases concern about climate change.  _Frontiers in Psychology, Environmental Psychology_ , 2019.
  * IMD [2018] IMD.  Report on the exceptionally heavy rainfall over kerala during 1st – 19th august 2018.  _Report on the exceptionally heavy rainfall over Kerala during 1st – 19th August 2018, Earth System science Organisation, Ministry of Earth Sciences_ , 2018.
  * Mishra and Shah [2019] Vimal Mishra and Harsh L. Shah.  Hydroclimatological perspective of the kerala flood of 2018.  _Journal Geological society of India, Vol.92_ , pages 645–650, 2019\.
  * Lal et al. [2020] Preet Lal, Aniket Prakash, Amit Kumar, Prashant K. Srivastava, Purabi Saikia, A.C. Pandey, Parul Srivastava, and M.L. Khan.  Evaluating the 2018 extreme flood hazard events in kerala, india.  _Remote sensing Letters, Vol. 11, No. 5_ , page 436–445, 2020.
  * Sudheer et al. [2019] K. P. Sudheer, S. Murty Bhallamudi, Balaji Narasimhan, Jobin Thomas, V. M. Bindhu, Vamsikrishna Vema, and Cicily Kurian.  Role of dams on the floods of august 2018 in periyar river basin, kerala.  _Current Science, Vol.116, No.5_ , 2019.
  * Hunt and Menon [2020] Kieran M. R. Hunt and Arathy Menon.  The 2018 kerala floods: a climate change perspective.  _Climate Dynamics, Vol. 54_ , page 2433–2446, 2020.
  * Kiran [2017] S. R. Kiran.  General circulation and principal wave modes in andaman sea from observations.  _Indian Journal of Science and Technology, Vol 10(24)_ , 2017.  doi:10.17485/ijst/2017/v10i24/115764.
  * Riehl [1954] H. Riehl.  _Tropical Meteorology; McGraw Hill, New York_ , 1954.
  * Charney [1969] J. G. Charney.  The intertropical convergence zone and the hadley circulation of the atmosphere.  _Proc. WMO/IUCG Symp. Numer. Weather Predict, Jpn. Meteorol. Agency_ , pages 73–79, 1969.
  * Gadgil [2018] Sulochana Gadgil.  The monsoon system: Land–sea breeze or the itcz?  _J. Earth Syst. Sci. (2018) 127:1,_ , pages 1–29, 2018.
  * Byrne et al. [2018] Michael P. Byrne, Angeline G. Pendergrass, Anita D. Rapp, and Kyle R. Wodzicki.  Response of the intertropical convergence zone to climate change:location, width, and strength.  _Current Climate Change Reports (2018), Vol. 4_ , page 355–370, 2018\.
  * Byrne and Thomasi [2019] Michael P. Byrne and Rhidian Thomasi.  Dynamics of itcz width: Ekman processes, non-ekman processes, and links to sea surface temperature.  _Journal of the Atmospheric Science, Vol. 76_ , page 2869–370, 2019\.
  * Murtugudde et al. [2007] Raghu Murtugudde, Richard Seager, and Prasad Thoppil.  Arabian sea response to monsoon variation.  _Paleoceanography, Vol. 22, PA4217_ , 2007.
  * Reed and Recker [1971] R. J. Reed and E. E. Recker.  Structure and properties of synoptic-scale wave disturbances in the equatorial western pacific.  _Journal of the Atmospheric Science, Vol. 28_ , 1971.
  * Ferreira and Schubert [1997] R. N. Ferreira and W. H. Schubert.  Barotropic aspects of itcz breakdown.  _Journal of the Atmospheric Science, Vol. 54_ , 1997.
  * Serra et al. [2008] Yolande L. Serra, George N. Kiladis, and Meghan F. Cronin.  Horizontal and vertical structure of easterly waves in the pacific itcz.  _Journal of the Atmospheric Science, Vol. 65_ , 2008.
  * Chen and Huang [2009] Guanghua Chen and Ronghui Huang.  Interannual variations in mixed rossby–gravity waves and their impacts on tropical cyclogenesis over the western north pacific.  _Journal of Climate, Vol. 22_ , 2009.
  * Raupp and Dias [2005] Carlos F. M. Raupp and Pedro L. Silva Dias.  Excitation mechanism of mixed rossby–gravity waves in the equatorial atmosphere: Role of the nonlinear interactions among equatorial waves.  _Journal of the Atmospheric Science, Vol. 62_ , 2005.
  * Dickinson and Molinari [2002] Michael Dickinson and John Molinari.  Mixed rossby–gravity waves and western pacific tropical cyclogenesis. part i: Synoptic evolution.  _Journal of the Atmospheric Science, Vol. 59_ , 2002.
  * Yang and Ingersol [2011] Da Yang and Andrew P. Ingersol.  Testing the hypothesis that the mjo is a mixed rossby–gravity wave packet.  _Journal of the Atmospheric Science, Vol. 68_ , 2011.
  * Wyrtki [1973] K. Wyrtki.  An equatorial jet in the indian ocean.  _Science, 181_ , pages 262–264, 1973.
  * Kiladis et al. [2009] G. N. Kiladis, M. C. Wheeler, P. T. Haertel, K. H. Straub, and P. E. Roundy.  Convectively coupled equatorial waves.  _Reviews of Geophysics, 47, RG2003_ , 2009.  doi:10.1029/2008RG000266.

Table 1: Rainfall at Pavannur & Mankadavu in Kerala during the period of anomalous rainfall in 2018 Period of Anomalous Rainfall | Pavannur Rainfall (mm/day) | Mankadavu Rainfall (mm/day)
---|---|---
Notation | Period | Max. | Anomaly | Max. | Anomaly
T1 | 08-June-2018 to 11-June-2018 | 98.50 | 47.41 | 69.20 | 51.94
| (Peak Rainfall on 09-June-2018) |  |  |  |
T2 | 03-July-2018 to 07-July-2018 | 105.70 | 90.06 | 7.27 | 4.30
| (Peak Rainfall on 05-July-2018) |  |  |  |
T3 | 16-July-2018 to 26-July-2018 | 196.50 | 148.80 | 49.52 | 11.01
| (Peak Rainfall on 18-July-2018) |  |  |  |
T4 | 04-August-2018 to 12-August-2018 | 249.50 | 149.30 | 97.26 | 54.81
| (Peak Rainfall on 08-August-2018) |  |  |  |
Figure 1: (top) TRMM daily accumulated precipitation (combined Microwave-IR)
estimate with gauge calibration performed over land (in millimetres) at
Mankadavu (Idukki District) and Pavannur (Kannur district) of Kerala in 2018.
(bottom) Daily rainfall anomalies (in millimetres) of 2018, measured with
respect to 7 year mean at Mankadavu and Pavannur. Figure 2: Spatial
distribution of TRMM daily accumulated precipitation (in millimetres) overlaid
with NCEP/NCAR Reanalysis surface winds (at 1000 hPa) on the days of peak
rainfall during T1, T2, T3 and T4. Figure 3: (a) The daily variation of
latitudinal position of Intertropical Convergence Zone (ITCZ) along $76^{o}$E.
(b) The meridional daily rainfall maximum inside ITCZ in millimeters along
$76^{o}$E. (c) The daily variation in the width of ITCZ in kilometres. Red
curve indicates the values corresponding to the year 2018 and blue curve
indicates the values averaged over 7 years. Figure 4: Variation of NCEP/NCAR
Reanalysis Zonal Wind Anomalies (in metres per second) with altitude (as
pressure levels in hecto pascals) over Kerala (averaged over $75^{o}$E to
$77^{o}$E & $8^{o}$N to $12^{o}$N). Figure 5: The Hovemöller diagram showing
temporal evolution of NCEP/NCAR Reanalysis Zonal Wind Anomalies (in metres per
second) at 700hPa pressure level along $10^{o}$N. The segment of the diagram
from $50^{o}$E to $150^{o}$E is enlarged for detailed inspection. Grey colored
arrows indicate phase propagation of signals in the westward direction. Figure
6: The Hovemöller diagram showing temporal evolution of NCEP/NCAR Reanalysis
Meridional Wind Anomalies (in metres per second) at 700hPa pressure level
along $10^{o}$N. The segment of the diagram from $50^{o}$E to $150^{o}$E is
enlarged for detailed inspection. Grey colored arrows indicate phase
propagation of signals in the westward direction and black colored arrow
indicates the eastward propagation of the wave group. $C_{p}$ and $C_{g}$
represents the phase and group velocities of waves respectively. Figure 7: The
variation of NCEP/NCAR Reanalysis Omega (indicative of Vertical velocity)
Anomalies (in Pa/s) over Kerala (averaged over $75^{o}$E to $77^{o}$E &
$8^{o}$N to $12^{o}$N) in 2018. Negative (positive) values represent rising
(sinking) of air parcels. Grey colored arrows indicate upward propagation of
the phase of waves and brown colored arrows indicate the propagation of wave
group. Figure 8: Spatial variation of NCEP/NCAR Reanalysis Relative vorticity
Anomalies (in $s^{-1}$) over the Equatorial Indian Ocean at 700 hPa level from
$16^{th}$ July to $20^{th}$ July 2018. The dashed red (blue) arrows indicate
the westward propagation of cyclonic (anti-cyclonic) vorticity anomalies.
Figure 9: Amplitude spectrum of NCEP/NCAR Reanalysis Relative vorticity
Anomalies (in $s^{-1}$) over equatorial Indian Ocean region (averaged over
$60^{o}$E to $100^{o}$E & $5^{o}$N to $20^{o}$N) at 700hPa level in 2018.
Figure 10: (first & second figures) 20 to 40 days bandpass-filtered NCEP/NCAR
Reanalysis Zonal & Meridional winds (in m/s) at 700hPa level averaged over the
latitudes $5^{o}$N to $15^{o}$N in 2018. (third & fourth figures) 5 to 8 days
bandpass-filtered NCEP/NCAR Reanalysis Zonal & Meridional winds (in m/s) at
700hPa level averaged over the latitudes $5^{o}$N to $15^{o}$N in 2018. Grey
colored arrows indicate propagation of phase in the westward direction.
$C_{p}$ represents the phase speed of waves. Figure 11: (top) Spatial
distribution of 20 to 40 days bandpass-filtered NCEP/NCAR Reanalysis Relative
vorticity Anomalies (in $s^{-1}$) at 700hPa level on the days of peak rainfall
during T1, T2, T3 and T4. (bottom) Spatial distribution of 5 to 8 days
bandpass-filtered NCEP/NCAR Reanalysis Relative vorticity Anomalies (in
$s^{-1}$) at 700hPa level for the same days. Positive (negative) values
indicate cyclonic (anti-cyclonic) vorticity in the Northern Hemisphere. Figure
12: (top) Spatial distribution of 20 to 40 days bandpass-filtered NCEP/NCAR
Reanalysis Omega Anomalies (in Pa/s) at 700hPa level on the days of peak
rainfall during T1, T2, T3 and T4. (bottom) Spatial distribution of 5 to 8
days bandpass-filtered NCEP/NCAR Reanalysis Omega Anomalies (in Pa/s) at
700hPa level for the same days. Negative (positive) values represent rising
(sinking) of air parcels.


# Gravitational Wave Cosmology: Be Careful of the Black Hole Mass Spectrum

Grégoire Pierra  0000-0003-3970-7970 Université Lyon, Université Claude
Bernard Lyon 1, CNRS, IP2I Lyon/IN2P3,UMR 5822, F-69622 Villeurbanne, France

###### Abstract

Gravitational waves (GWs) from compact binary coalescences (CBCs) offer
insights into the universe’s expansion. The “spectral siren” method, used
without electromagnetic counterparts (EMC), infers cosmic expansion (Hubble
constant $H_{0}$) by relating detector and source frame masses of black hole
(BH) mergers. However, heuristic mass models (broken power law, power law plus
peak, multipeak) may introduce biases in $H_{0}$ estimation, potentially up to
3$\sigma$ with 2000 detected GW mergers. These biases stem from models’
inability to consider redshift evolution and unexpected mass features. Future
GW cosmology studies should employ adaptable source mass models to address
these issues.

##  1 Introduction

Since September 2015, the LIGO-Virgo-KAGRA collaboration has detected over 90
GW signals from CBCs [1, 2, 3], offering direct measurements of luminosity
distances ($d_{L}$) without relying on specific cosmological models. However,
without an EMC, determining the redshift of the source to measure the
universe’s expansion rate becomes a complex task, a topic currently debated
[4]. To address this, methods like the spectral siren technique have emerged,
estimating redshift indirectly using the intrinsic relationship between
redshift and masses in the detector and source frames [5]. The accuracy of
this approach relies on chosen phenomenological models for BBH source mass
distribution. This study investigates how well common source mass models fit
complex BBH distributions, aiming to uncover potential biases in $H_{0}$.
Using synthetically generated BBH populations and hierarchical Bayesian
inference [6], we simultaneously measure $H_{0}$ alongside phenomenological
BBH population mass model.

##  2 Method

###  2.1 Simulation of GW observations

The workflow for generating GW observations involves calculating an
approximate optimal signal-to-noise ratio (SNR) from a BBH merger population,
characterized by source masses and redshift. Instead of using the full
waveform and matched filter, we use a proxy formula for the optimal SNR [7,
8]:

$\rho=\delta\times\rho_{9}\left[\frac{\mathcal{M}_{c}}{\mathcal{M}_{c,{\rm
9}}}\right]^{\frac{5}{6}}\left[\frac{d_{L,{\rm 9}}}{d_{L}}\right],$ (1)

where $\mathcal{M}_{c}$ is the chirp mass, $d_{L}$ is the luminosity distance,
$\rho$ is the optimal SNR, and $\delta$ is a projection factor. We assume a
flat $\Lambda CDM$ cosmology with $H_{0}=67.8$ ${\rm km\,s^{-1}\,Mpc^{-1}}$
and $\Omega_{m,0}=0.308$, based on Planck measurements [9]. The reference
chirp mass, luminosity distance and optimal SNR are set to
$\mathcal{M}_{c,{\rm 9}}=25$ $\rm{M}_{\odot}$, $d_{L,{\rm 9}}=1.5$ $\rm Gpc$
and $\rho_{9}=0$ [10]. The projection factor $\delta$ accounts for the
3-detector network’s sensitivity, including GW polarization and sky
localization. To simulate detector noise, we calculate a ”detected” SNR
$\rho^{\rm det}$ from the optimal SNR using a $\chi^{2}$ distribution. A GW is
considered detected if $\rho^{\rm det}$ exceeds 12 and the GW frequency from
the innermost stable circular orbit (ISCO) is above 15 Hz. For this study, we
assume perfect measurement of detector frame masses and luminosity distances
to reduce computational load and focus on systematic effects from
phenomenological mass models.

###  2.2 Bayesian inference

The detection of GW events is modeled as an inhomogeneous Poisson process with
selection biases [11, 6, 12, 8]. For $N_{\rm GW}$ detected GW signals over an
observation time $T_{\rm obs}$, the probability of obtaining a specific GW
dataset {x} given population hyperparameters $\Lambda$ is:

$\displaystyle\mathcal{L}(\\{x\\}|\Lambda)$ $\displaystyle\propto$
$\displaystyle e^{-N_{\rm exp}(\Lambda)}\prod_{i}^{N_{\rm GW}}T_{\rm
obs}\int{\rm d}\theta{\rm d}z\mathcal{L}_{\rm
GW}(x_{i}|\theta,z,\Lambda)\frac{1}{1+z}\frac{{\rm d}\rm N_{\rm CBC}}{{\rm
d}\theta{\rm d}z{\rm d}t_{s}}(\Lambda).$ (2)

The hyperparameters $\Lambda$ include BBH population parameters and
cosmological expansion parameters (e.g., $H_{0}$ and $\Omega_{\rm m}$). The
individual GW likelihood $\mathcal{L}_{\rm GW}(x_{i}|\theta,\Lambda)$ accounts
for errors in intrinsic parameters $\theta$, specifically source frame masses.
The CBC rate, parameterized for the spectral siren method, is:

$\frac{{\rm d}\rm N_{\rm CBC}}{{\rm d}\theta{\rm d}z{\rm
d}t_{s}}=R_{0}\psi(z;\Lambda)p_{\rm pop}(\vec{m}_{\rm
s}|\Lambda)\times\frac{{\rm d}V_{c}}{{\rm d}z},$ (3)

where $R_{0}$ is the local merger rate, $\vec{m}_{\rm s}=(\rm m^{\rm
s}_{1},m^{\rm s}_{2})$ are the source frame masses, and $\psi(z;\Lambda)$ is
the merger rate evolution with redshift. The term $\rm p_{\rm
pop}(\vec{m}_{s}|\Lambda)$ is the probability density function for the source
frame masses and $\frac{{\rm d}V_{c}}{{\rm d}z}$ is the differential of the
comoving volume. We assume the CBC mass spectrum does not evolve with
redshift, a common assumption in the literature [13, 14, 6, 12]. The
implications of this assumption on $H_{0}$ bias are discussed subsequently.

###  2.3 Population models of BBHs

Phenomenological models for CBC merger rates, masses, and redshift play a
pivotal role in this study. Current spectral siren analyses, exemplified by
[13], typically employ simple parametric models for the distribution of source
frame masses and CBC merger rates. We utilize three mass models: Broken Power
Law (BPL), Power Law plus Peak (PLP), and Multipeak (MLTP). The BPL [15] is
the most straightforward, featuring a power law with low mass smoothing and a
break at $\rm m_{\rm break}$, accounting for effects of stellar progenitor
metallicity and the pair-instability supernovae gap. The PLP [16] expands on
the BPL by adding a Gaussian peak to capture a potential BBH accumulation
before the pair instability supernovae gap (PISN), with parameters like peak
position and width allowing for flexibility. The MLTP model further extends
the PLP by introducing a second Gaussian peak at higher masses, accommodating
BBHs from second-generation black holes due to hierarchical mergers.

##  3 Results

###  3.1 Sanity checks

We simulate three test populations of BBH mergers, each with different source
mass distributions and CBC merger rate. For each population, we generate 2000
detected GW events for use in spectral siren analysis. This number of
detections typically ensures a data-informed posterior within the prior range
$H_{0}\in[20,140]{\rm km\,s^{-1}\,Mpc^{-1}}$. We employ hierarchical Bayesian
inference to jointly estimate the cosmological parameter $H_{0}$ alongside all
mass and CBC merger rate parameters. Left plot of Fig. 1 presents marginal
posterior distributions for $H_{0}$ obtained for each source mass model. The
true $H_{0}$ value is estimated within the 90% C.I. for all simulations, with
true values of other population parameters also contained in their respective
90% C.I. The PLP and MLTP mass models yield more precise $H_{0}$ constraints
than the BPL, owing to sharper features in the source mass spectrum.

Figure 1: Left: Marginal posteriors of $H_{0}$ obtained via spectral sirens
inference. The red line is the true value injected. Right:Reconstructed
spectrum of the primary mass in source frame from the three inferences.

Right plot of Fig. 1 displays posterior predictive checks (PPCs) for the BPL,
PLP, and MLTP mass models. These PPCs demonstrate accurate reconstruction of
mass spectrum structures by the inference, including all features present in
the mass spectra. Thus, when using the correct mass model for inference, both
population and cosmological expansion parameters are accurately recovered.
Additionally, this test demonstrates the importance of features in the BBH
mass spectrum for improving precision in cosmological expansion parameters.

###  3.2 Linear redshift evolution

We conduct a test to assess the robustness of redshift-independent mass models
to a potential evolution of the mass spectrum with redshift. We assume a
linear dependence between mass spectrum features and redshift. Here, we modify
the PLP mass model to include a linearly evolving peak position with redshift,
while keeping the edges fixed. This evolution could arise from interplay
between PISN mass scale and binary formation-to-merger time-delays [17, 18].
The modified PLP model incorporates redshift evolution as follows:

$\mu_{g}(z)=\mu_{g}^{0}+z\left(\mu_{g}^{1}-\mu_{g}^{0}\right),$ (4)

where $\mu_{g}^{0}$ and $\mu_{g}^{1}$ denote peak positions at $z=0$ and $z=1$
respectively, with $\mu_{g}^{0}=30\rm{M}_{\odot}$ and varying $\mu_{g}^{1}$
from $25\rm{M}_{\odot}$ to $35\rm{M}_{\odot}$. Each simulation involves
hierarchical Bayesian inference using the redshift-independent PLP to
reconstruct $H_{0}$ and other population parameters from 2000 GW detections.
Left of Fig. 2 shows $H_{0}$ estimations with 68.3% C.I. for these analyses,
indicating $H_{0}$ variations as $\mu_{g}^{1}$ changes. Even a mild mass
spectrum evolution biases the inferred $H_{0}$ significantly, linearly
proportional to the redshift evolution magnitude. This bias arises from
inaccurate reconstruction of the source-frame mass spectrum, consistent across
independent simulation realizations. Right of Fig. 2 presents PPC for this
test. The PLP model fails to capture this evolution, consistently
reconstructing a peak at 30 $\rm{M}_{\odot}$. Misplacement of the Gaussian
component at higher redshifts induces $H_{0}$ bias: for $\mu_{\rm
g}^{1}>\mu_{\rm g}^{0}$, $H_{0}$ is biased higher, and vice versa. This
suggests potential bias in GW-based $H_{0}$ estimation from redshift-evolving
mass spectrum features.

Figure 2: Left: Inferred value of $H_{0}$ as a function of the position of the
Gaussian peak at $z=1$. Right:Reconstructed spectrum of the primary mass in
source frame from the three inferences.

###  3.3 Complex BBHs population

Utilizing an astrophysically motivated synthetic BBH catalog (A03) [19], we
generate GW catalogs containing $2000$ detected events. Each catalog undergoes
full hierarchical Bayesian inference to estimate posteriors on $H_{0}$ and
other population parameters, using redshift rate and three redshift-
independent BBH mass models. Left of Fig. 3 depicts the marginal posterior
distributions of estimated Hubble constants. Both PLP and MLTP mass models
exhibit an $H_{0}$ bias towards higher values. Specifically, PLP yields
$H_{0}=103^{+11}_{-10}{\rm km\,s^{-1}\,Mpc^{-1}}$ and $H_{0}=87^{+9}_{-8}{\rm
km\,s^{-1}\,Mpc^{-1}}$ at 68.3% C.I., excluding the true value of $H_{0}$ from
the 99.7% and 95% C.I. respectively. In contrast, the BPL mass model yields an
uninformative $H_{0}$ posterior, encompassing the injected $H_{0}$ within the
68% C.I. The lack of strong mass features in BPL contributes to this less
informative posterior. PLP and MLTP underestimate the presence of BBHs between
$40\rm{M}_{\odot}$ and $80\rm{M}_{\odot}$, leading to an overestimation of
$H_{0}$ by pushing events to higher redshifts with lower source masses.
Although PLP and MLTP possess enough flexibility to approximate smooth mass
distributions akin to the A03 catalog, they fail to accurately reconstruct the
mass spectrum and consequently introduce $H_{0}$ bias.

Figure 3: Left: Marginal posterior of $H_{0}$ with the A03 BBH population.
Right:Marginal posterior of $H_{0}$ with the A03 BBH population, where the
redshift evolution has been blinded to the masses.

To further explore the potential impact of redshift evolution on the mass
spectrum and its influence on $H_{0}$ bias, we conducted a simulation by
blinding the A03 mass spectra to redshift evolution. This involved randomly
shuffling pairs of BBH merger redshifts and masses, removing the redshift
dependency while preserving the mass spectrum’s shape. We focused on the PLP
and MLTP mass models due to their notable $H_{0}$ bias. Right of Fig. 3
presents the joint inference results with the PLP mass model, yielding
$H_{0}=67.0^{+9.7}_{-8.9}{\rm km\,s^{-1}\,Mpc^{-1}}$. Similar results were
obtained with MLTP ($H_{0}=75.13^{+7.8}_{-7.9}{\rm km\,s^{-1}\,Mpc^{-1}}$),
with the true $H_{0}$ value falling within the 68.3% C.I. of the posterior.
Our test confirms that redshift evolution of the mass spectrum is the
significant contributor to $H_{0}$ bias.

##  4 Conclusion

This study delves into systematic biases affecting $H_{0}$ estimation and BBH
mass spectrum reconstruction in spectral siren Cosmology. Through simulations
with various mass models and synthetic BBH catalogs, we uncovered critical
insights. We revealed that simplistic models lacking sharp features may bias
$H_{0}$ if true BBH populations exhibit local over-densities. Additionally,
neglecting redshift evolution in source frame mass features can skew $H_{0}$
estimates, as underestimated mass spectra shift $H_{0}$ higher, and vice
versa. Using the A03 BBH catalog, we found $H_{0}$ biases of up to $3\sigma$
for PLP and MLTP inferences. Removing redshift evolution eliminated biases.
Future GW cosmology studies should employ flexible mass models capable of
accommodating such complexities to avoid biased $H_{0}$ estimations,
especially with numerous detected GW events.

## Acknowledgments

The authors are grateful for computational resources provided by the LIGO
Laboratory and supported by the National Science Foundation Grants PHY-0757058
and PHY-0823459. We thank Viola Sordini for the helpful comments and
discussions on the manuscript. MM acknowledges financial support from the
European Research Council for the ERC Consolidator grant DEMOBLACK, under
contract no. 770017, and from the German Excellence Strategy via the
Heidelberg Cluster of Excellence (EXC 2181 - 390900948) STRUCTURES.

## References

  * [1] R. Abbott et al.  GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo during the Second Part of the Third Observing Run.  Phys. Rev. X, 13(4):041039, 2023.
  * [2] R. Abbott et al.  GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run.  Phys. Rev. X, 11:021053, 2021.
  * [3] R. Abbott et al.  GWTC-2.1: Deep Extended Catalog of Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run.  arXiv, 8 2021.
  * [4] Wendy L. Freedman.  Measurements of the Hubble Constant: Tensions in Perspective.  Astrophys. J., 919(1):16, 2021.
  * [5] Jose María Ezquiaga and Daniel E. Holz.  Spectral Sirens: Cosmology from the Full Mass Distribution of Compact Binaries.  Phys. Rev. Lett., 129(6):061102, 2022.
  * [6] Simone Mastrogiovanni, Grégoire Pierra, Stéphane Perriès, Danny Laghi, Giada Caneva Santoro, Archisman Ghosh, Rachel Gray, Christos Karathanasis, and Konstantin Leyde.  ICAROGW: A python package for inference of astrophysical population properties of noisy, heterogeneous and incomplete observations.  arXiv, 5 2023.
  * [7] Maya Fishbach, Daniel E. Holz, and Will M. Farr.  Does the Black Hole Merger Rate Evolve with Redshift?  Astrophys. J. Lett., 863(2):L41, 2018.
  * [8] S. Mastrogiovanni, K. Leyde, C. Karathanasis, E. Chassande-Mottin, D. A. Steer, J. Gair, A. Ghosh, R. Gray, S. Mukherjee, and S. Rinaldi.  On the importance of source population models for gravitational-wave cosmology.  Phys. Rev. D, 104(6):062009, 2021.
  * [9] P. A. R. Ade et al.  Planck 2015 results. XIII. Cosmological parameters.  Astron. Astrophys., 594:A13, 2016.
  * [10] B. P. Abbott et al.  Prospects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo and KAGRA.  Living Rev. Rel., 21(1):3, 2018.
  * [11] Salvatore Vitale, Davide Gerosa, Will M. Farr, and Stephen R. Taylor.  Inferring the properties of a population of compact binaries in presence of selection effects.  arXiv, 7 2020.
  * [12] Simone Mastrogiovanni, Danny Laghi, Rachel Gray, Giada Caneva Santoro, Archisman Ghosh, Christos Karathanasis, Konstantin Leyde, Daniele A. Steer, Stephane Perries, and Gregoire Pierra.  A novel approach to infer population and cosmological properties with gravitational waves standard sirens and galaxy surveys.  arXiv, 5 2023.
  * [13] R. Abbott et al.  Constraints on the Cosmic Expansion History from GWTC–3.  Astrophys. J., 949(2):76, 2023.
  * [14] Rachel Gray et al.  Joint cosmological and gravitational-wave population inference using dark sirens and galaxy catalogues.  arXiv, 8 2023.
  * [15] R. Abbott et al.  Population Properties of Compact Objects from the Second LIGO-Virgo Gravitational-Wave Transient Catalog.  Astrophys. J. Lett., 913(1):L7, 2021.
  * [16] Colm Talbot and Eric Thrane.  Measuring the binary black hole mass spectrum with an astrophysically motivated parameterization.  Astrophys. J., 856(2):173, 2018.
  * [17] Suvodip Mukherjee.  The redshift dependence of black hole mass distribution: is it reliable for standard sirens cosmology?  Mon. Not. Roy. Astron. Soc., 515(4):5495–5505, 2022.
  * [18] Christos Karathanasis, Benoît Revenu, Suvodip Mukherjee, and Federico Stachurski.  GWSim: Python package for creating mock GW samples for different astrophysical populations and cosmological models of binary black holes.  A&A, 677:A124, September 2023.
  * [19] Michela Mapelli, Yann Bouffanais, Filippo Santoliquido, Manuel Arca Sedda, and M. Celeste Artale.  The cosmic evolution of binary black holes in young, globular, and nuclear star clusters: rates, masses, spins, and mixing fractions.  Mon. Not. Roy. Astron. Soc., 511(4):5797–5816, 2022.