# Homogenization and trend estimation of GNSS water vapour observations and atmospheric reanalyses. Application to climate change analysis.

Khanh-Ninh Par

► To cite this version:

Khanh-Ninh Par. Homogenization and trend estimation of GNSS water vapour observations and atmospheric reanalyses. Application to climate change analysis. Environmental Sciences. IPGP, Université Paris Cité, 2023. English. NNT: tel-04354790

HAL Id: tel-04354790

https://theses.hal.science/tel-04354790

Submitted on 19 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
Université Paris Cité

Institut de Physique du Globe de Paris

École doctorale Sciences de la Terre et de l'environnement et physique de l'Univers

STEP'UP n°560

Équipe de géodésie

# Homogénéisation et estimation des tendances des observations GNSS de la vapeur d'eau et des réanalyses atmosphériques. Application à l'analyse du changement climatique.

Homogenization and trend estimation of GNSS water vapour observations and atmospheric
reanalyses. Application to climate change analysis.

Par Khanh-Ninh Nguyen

Thèse de doctorat de Sciences de la Terre

Dirigée par Olivier Bock
et Emilie Lebarbier

Présentée et soutenue publiquement le 8 Décembre 2023

Devant un jury composé de :

Mme Hélène Brogniez, PU, LATMOS, Université de Versailles Saint-Quentin-en-Yvelines, rapporteur

M. Athanassios Argiriou, Full professor, University of Patras, rapporteur

M. Hervé Douville, Directeur de Recherche, Météo-France, Université Paul Sabatier, Examinateur

M. Clément Narteau, Professeur, IPGP, Université Paris Cité, Examinateur

Mme Emilie Lebarbier, PU, Université Paris Nanterre, co-directrice de thèse

M. Olivier Bock, Directeur de Recherche, IGN, IPGP, Université Paris Cité, directeur de thèse
<figure><img src="image_1.png" /></figure>
Acknowledgment

First and foremost, I extend my deepest gratitude to my supervisors, Olivier Bock, Emilie Lebarbier, and Françoise Guichard, for their invaluable guidance, patience, and unwavering support throughout my PhD journey, from preparation to completion. Special appreciation goes to Olivier, a professional researcher, for his consistent and continued support and tolerance. His countless feedback and meticulous guidance have been pivotal in refining my academic writing and sharpening my arguments.

I am deeply thankful to Emilie Lebarbier for her rigorous mathematical perspective and her wealth of creative ideas. Our weekly meetings kept me on track and progressively shaped my academic growth. Francoise Guichard, though no longer with us, provided warm support at the outset of my PhD journey during the challenging times of COVID quarantine. Her guidance in both scientific and French cultural aspects has left an indelible mark on my development.

I also could not have undertaken this journey without my PhD defense committee members: Prof. Athanassios Argiriou, Prof. Hélène Brogniez, Dr. Hervé Douville, and Prof. Clément Narteau. I am grateful for your acceptance of the invitation to be part of my committee, and for your insightful feedback and constructive critiques, which have significantly enhanced this thesis.

My appreciation goes out to the Institut de Physique du Globe de Paris (IPGP) and Institut national de l'information géographique et forestière (IGN) for their financial support and for providing an environment conducive to research and innovation. This opportunity has been pivotal in my academic and personal development.

Additionally, I owe a great deal to my wonderful lab mates, whose support has been a constant source of motivation. A heartfelt thanks to Maylis, Ana, and Marie for welcoming me warmly to the team, sharing your experiences, and providing encouragement and support both academically and personally. My gratitude also extends to Louis-Marie, Charlotte, Laetitia, Boris, Jiao for their companionship throughout this journey, and to Paul, Kevin, Julien, Kristel, Laurent, and Gwendoline for their invaluable support and advice in research and academic matters. Thanks also to Anne and Nora for their assistance with administrative tasks.

Lastly, but most importantly, a heartfelt thank you to my family. Your unwavering love, understanding, and encouragement have been the pillars of my strength. To my mom, dad, Duc, and Quan: Thank you for everything. This PhD thesis is dedicated to you.
Résumé (en français)

La vapeur d'eau joue un rôle important dans le bilan énergétique et le cycle hydrologique de la Terre. Avec le réchauffement climatique, la vapeur d'eau dans l'atmosphère a tendance à augmenter, mais la répartition géographique des tendances de la vapeur d'eau dans les modèles climatiques et dans les observations reste assez incertaine. Les estimations de la vapeur d'eau intégrée (IWV) fournies par le système mondial de navigation par satellite (GNSS) constituent une nouvelle source de données qui peut servir à confirmer les tendances observées et simulées de la vapeur d'eau. Des études précédentes ont montré l'existence d'inhomogénéités dans les enregistrements IWV du GNSS dues aux changements des instruments de mesures des stations et dans les méthodes de traitement des données, et ont recommandé leur homogénéisation pour les applications climatiques.

Comme les approches récentes d'homogénéisation de données climatiques, cette étude utilise une méthode statistique de segmentation pour détecter les changements abrupts appelées ruptures dans la moyenne des séries de différences IWV entre les données GNSS et les données de réanalyse. Les séries de différences sont utilisées pour atténuer la forte variabilité temporelle inhérente aux séries IWV, rendant difficile la détection de petits sauts sur les séries brutes. Cette étude utilise le package R GNSSseg développé par Quarello (2020).

Le premier objectif de cette thèse était d'étudier la sensibilité des résultats de la segmentation et des estimations de tendances aux propriétés des données d'entrée. Deux jeux de données GNSS spécifiques ont été considérés (IGS repro1 et CODE REPRO2015) ainsi que deux réanalyses différentes (ERA-Interim et ERA5). Il s'avère que seuls 45 à 50% de ruptures communes sont détectées entre les deux jeux de données GNSS. La longueur de la série ou l'utilisation de données auxiliaires dans le traitement des données GNSS IWV a un impact plus faible, avec 70 – 80% de détections similaires. Les résultats de la segmentation sont sensibles au niveau de bruit et à la présence d'un biais périodique dans les données, qui sont principalement dus à des différences de représentativité entre les observations GNSS et les réanalyses. Les données récentes de GNSS et de réanalyses présentent des différences moins importantes et permettent la détection de sauts plus faibles. Les estimations des tendances résultantes sont sensibles au nombre et à la position des ruptures. C'est pourquoi la validation des ruptures détectées est une étape cruciale avant l'estimation des tendances. Il a été constaté que les changements d'instruments GNSS disponibles sous forme de métadonnées peuvent expliquer 35% des ruptures, en laissant 65% non documentés ou correspondants à des inhomogénéités dans les réanalyses.

Le deuxième objectif de cette thèse était de développer une méthode de classification automatique de post-détection pour distinguer les ruptures du GNSS de la réanalyse. La méthode "d'attribution" proposée combine les données GNSS et de réanalyse d'une station principale avec des données similaires de stations voisines. Chaque paire "station principale-station voisine" comprend quatre séries de base à partir desquelles six séries de différences sont formées. Ensuite, pour chaque rupture détectée dans la série principale, un test de significativité du saut de moyenne à cet instant est effectué, et une règle statistique de prédiction est construite pour attribuer la rupture testée au GNSS ou à la réanalyse. Lorsque plusieurs stations proches sont disponibles, une solution combinée est proposée. L'aspect original de la méthode d'attribution développée est l'utilisation d'une méthode d'inférence généralisée des moindres carrés, pour tenir compte de l'hétéroscédasticité et de l'autocorrélation dans les données, et d'une méthode d'apprentissage automatique. La principale nouveauté réside dans la manière dont la règle de prédiction est construite, en utilisant les résultats des tests obtenus à partir des données réelles
avec une procédure de rééchantillonnage. Une analyse de sensibilité montre que les prédictions dépendent d'une manière ou d'une autre de la stratégie de rééchantillonnage, en utilisant des échantillons équilibrés ou déséquilibrés, du niveau de significativité du test et également du niveau de bruit dans les données. Un facteur limitant dans les réseaux éparses est le bruit élevé quand les séries sont très éloignées. Lorsque la méthode est appliquée aux données GNSS CODE REPRO2015 et à la réanalyse ERA5, 62% des ruptures testées sont attribuées au GNSS, 19% à la réanalyse et 10% à des changements à la fois dans le GNSS et la réanalyse.

**Mots-clés** : vapeur d'eau intégrée, GNSS, homogénéisation, tendance, segmentation, attribution, réanalyse, régression linéaire généralisée, classification supervisée
**Abstract** (in english)

Water vapour plays an important role in the Earth's energy balance and hydrological cycle. As the climate warms, water vapour in the atmosphere tends to increase, but the geographical distribution of water vapour trends in climate model simulations and observations remains quite uncertain. Global Navigation Satellite System (GNSS) estimates of Integrated Water Vapor (IWV) are a new data source that may serve for the verification of observed and simulated water vapour trends. Previous studies have reported inhomogeneities in the GNSS IWV records due to changes in the station instrumentation and in data processing methods, and recommended their homogenization for climate applications.

Following modern approaches to climate data homogenization, this study uses a statistical segmentation method to detect abrupt changes, called change-points, in the mean (or jumps) of IWV differences series between GNSS and reanalysis data. Difference series are used to mitigate the strong temporal variability inherent in the IWV series which makes the detection of small jumps difficult on the raw series. This study employs the GNSSseg R package developed by Quarello (2020).

The first objective of this thesis was to investigate the sensitivity of the segmentation results and trend estimates to the input data properties. Two specific GNSS datasets were considered as inputs (IGS repro1 and CODE REPRO2015) and two different reanalysis products (ERA-Interim and ERA5). It is found that only 45-50% similar change-points are detected when the input data are altered. Altering the series length or the auxiliary data used in the processing of GNSS IWV data has a smaller impact, with 70-80% similar detections. Segmentation results are sensitive to the noise magnitude and the presence of a periodic bias in the input data, mainly due to representativeness differences between the GNSS observations and reanalysis. More recent GNSS and reanalysis products have smaller differences and allow the detection of smaller jumps. The subsequent trend estimates are sensitive to the number and position of change-points. Therefore, validation of detected change-points is a crucial step before trend estimation. It is found that changes in GNSS instrumentation available as metadata can explain 35% of change-points, leaving 65% as either undocumented GNSS changes or changes in the reanalysis data.

The second objective of this thesis was to develop an automatic classification method, operated as a post-processing step, to distinguish GNSS and reanalysis change-points. The proposed "attribution" method combines GNSS and reanalysis data from a main station with similar data from nearby stations. Each main-nearby pair comprises four base series from which six series of differences are formed. Then, for each detected change-point in the main series, a significance test is applied to the mean before and after the change-point, and a classifier or a statistical predictive rule is constructed to attribute the tested change-point to GNSS or reanalysis. When several nearby stations are available, a combined solution is proposed. The original aspect of the developed attribution method is the use of a generalized least-squares inference method, to account for heteroscedasticity and autocorrelation in the data, and a machine learning classifier. The main novelty is the way in which the predictive rule is constructed, using test results obtained from the real data with a resampling procedure. Sensitivity analysis shows that the prediction results depend somehow on the resampling strategy, using balanced or imbalanced samples, on the test significance level, and also on the magnitude of the noise in the data. A limiting factor in spare networks is the large noise with long distance between the main and nearby stations.
When applied to the CODE REPRO2015 GNSS data and ERA5 reanalysis, 62% of the tested change-points are attributed to GNSS, 19% to the reanalysis, and 10% to changes in both GNSS and reanalysis.

**Keywords:** integrated water vapor, GNSS, homogenization, trend, segmentation, attribution, reanalysis, generalized least square, supervised classification
# Résumé substantiel (en français)

## Introduction

La vapeur d'eau joue un rôle important dans le bilan énergétique et le cycle hydrologique de la Terre. En tant que gaz à effet de serre, elle amplifie le réchauffement de la planète par son mécanisme de rétroaction positive. La compréhension et la quantification de l'évolution de la vapeur d'eau dans le contexte du changement climatique représentent un défi scientifique majeur, comme le soulignent des recherches récentes (Colman and Soden, 2021; Douville et al., 2022). À mesure que le climat se réchauffe, la vapeur d'eau dans l'atmosphère a tendance à augmenter, mais la répartition géographique des tendances de la vapeur d'eau dans les simulations de modèles climatiques et les observations reste assez incertaine (Hersbach et al., 2020; Flato et al., 2013). Dans ce contexte, le système mondial de navigation par satellite (GNSS) offre des estimations de la vapeur d'eau intégrée (IWV) avec une haute résolution temporelle disponible dans toutes les conditions météorologiques (Guerova et al., 2016), qui peuvent servir de source indépendante de données pour la vérification des tendances de la vapeur d'eau observées et simulées (Parracho, 2017). Néanmoins, l'homogénéité des données GNSS est sensible à divers facteurs, y compris aux changements des instruments de mesures, des méthodes de traitement des données, et environnementaux (Vey et al., 2009; Bock et al., 2010; Ning et al., 2016). Dans ce contexte, cette thèse vise à contribuer à l'effort d'homogénéisation des enregistrements GNSS IWV, ce qui est essentiel pour permettre des estimations précises des tendances de la vapeur d'eau atmosphérique.

## Bilan et objectifs

L'homogénéisation fait référence au processus d'identification et de correction des inhomogénéités (artificielles) dans des séries de données. Suivant les approches modernes de l'homogénéisation des données climatiques, une méthode statistique de segmentation pour détecter les changements abrupts, appelées ruptures, dans la moyenne des séries de différences IWV entre les données GNSS et les données de réanalyse (approche qualifiée de "relative") a été utilisée. L'utilisation de séries de différences est stratégique pour atténuer la forte variabilité temporelle inhérente aux séries IWV, rendant difficile la détection de petits sauts sur les séries brutes. La méthode de segmentation utilisée dans cette étude a été développée par Quarello (2020) et dédiée à l'analyse de ces séries de différences (en s'adaptant à ces caractéristiques). Elle a été validée par des études de simulations et a été reconnue comme l'une des plus efficaces lorsqu'elle a été comparée à d'autres méthodes sur des données benchmark dans l'étude de Van Malderen et al. (2020). Toutefois, la robustesse, la sensibilité et les limites de cette méthode par rapport à diverses caractéristiques des données et ses effets sur l'estimation des tendances n'avaient pas été étudiées. Il s'agit d'un aspect crucial pour son application pratique à des ensembles de données réelles. Le premier objectif de cette thèse est donc de combler cette lacune. Le second est de déterminer les origines des ruptures, qui peuvent provenir de données GNSS ou de données de réanalyse. Cette étape est appelée attribution. L'attribution précise des ruptures à leur source est essentielle avant de corriger les données IWV brutes de leurs sauts. Représentant un défi fondamental de l'approche relative, ce problème d'attribution demeure une limite clé des méthodes d'homogénéisation existantes. En réponse, cette thèse introduit une méthode automatique en post-traitement à la segmentation. Une étude poussée à la fois sur des simulations mais aussi sur les vraies données est menée afin de bien comprendre les effets des différents facteurs mis en jeu lors de cette étape d'attribution. Enfin, le dernier objectif est de corriger les séries temporelle GNSS IWV des sauts associes aux ruptures qui leurs sont attribuées.
## Sensibilité de la méthode de segmentation aux propriétés des données et étude de l'impact de ces propriétés sur les estimations de tendances

Une importante étude de la sensibilité de la méthode de segmentation aux propriétés des données a été réalisée, ainsi que l'évaluation de l'impact de ces propriétés sur les estimations de tendances. Cette analyse a comparé des paires de jeux de données selon quatre facteurs critiques: les méthodes de traitement des données GNSS (IGSrepro1 vs. CODE REPRO2015), la longueur des séries (17 ans vs. 25 ans), les données auxiliaires utilisées dans la conversion de la vapeur d'eau intégrée (IWV), et les sources des données de référence (ERA-Interim et ERA5). L'accent a été mis sur des propriétés spécifiques des données, telles que les valeurs moyennes, les niveaux de bruit et les biais périodiques, car ces propriétés ont un impact sur les résultats de la segmentation, en particulier en termes d'influence sur le nombre et le positionnement des ruptures. Les résultats indiquent des impacts significatifs sur les résultats de la segmentation lorsque l'on modifie le traitement GNSS et la réanalyse de référence. Seules 45 à 49% des ruptures sont similaires dans ces cas, contre 71 à 81% pour les deux autres facteurs (la longueur de la série et les données auxiliaires). Les changements notables dans le traitement GNSS comprennent des changements dans la correction ZHD a priori, le modèle d'étalonnage de l'antenne/radome et la fonction de cartographie. Les améliorations apportées au traitement CODE permettent de réduire le bruit et le biais périodique. De même, le passage de l'ERA-Intérim à l'ERA5 en tant que référence réduit les erreurs de représentativité, ce qui entraîne là aussi une réduction du bruit et des biais périodiques, facilitant ainsi la détection des « petites » ruptures. Le taux de validation des ruptures détectées par les métadonnées est entre 30% et 35% pour tous les ensembles de données, ce qui signifie que l'impact des propriétés des données sur ce taux est faible. Ce résultat suggère 65% des ruptures détectées sont soit des changements GNSS non documentés, soit des changements dans les données de réanalyse.

L'impact sur les estimations des tendances est étudié selon deux aspects. Premièrement, l'influence des quatre facteurs précédents sur les estimations des tendances de l'IWV est étudiée. Il apparaît que la modification de la longueur de la série temporelle a l'impact le plus fort, à la fois sur la moyenne et sur la dispersion des estimations des tendances à travers le réseau. Le changement de moyenne semble indiquer l'intensification du cycle de l'eau au cours de la décennie 2010-2020 par rapport à la décennie précédente. La dispersion réduite de la tendance avec la période plus longue est principalement le résultat de l'erreur standard plus petite des estimations. Deuxièmement, l'effet de la correction des inhomogénéités sur les estimations des tendances est étudié. Cet effet est particulièrement marqué lorsque seules les ruptures validées par les métadonnées GNSS (changements de récepteur, d'antenne, de radôme) sont corrigées. Cet impact est évident dans les tendances moyennes globales, la dispersion et les différences RMS par rapport à ERA5. La correction utilisant toutes les ruptures a également un impact sur les estimations des tendances, mais dans une mesure légèrement moindre. Ceci met en exergue la nécessité d'une étape d'attribution afin de retenir uniquement les ruptures attribuables au GNSS pour la correction. Enfin, la dispersion associée aux différentes estimations homogénéisées des tendances peut être utilisée comme mesure de l'incertitude dans l'estimation des tendances à une seule station. Elle s'élève à 0.1 – 0.2 kg m⁻² décennie⁻¹, ou 0.5-1% décennie⁻¹, ce qui confirme la faisabilité de la détection des tendances climatiques mondiales et régionales pertinentes avec les données GNSS IWV.
Développement de la méthode d'attribution

Le second objectif de la thèse porte sur le développement d'une méthode automatique d'attribution des ruptures. Elle a pour objectif de déterminer l'origine des ruptures détectées par la segmentation, i.e. de distinguer les ruptures du GNSS de celles de la réanalyse. Cette méthode consiste à combiner les données GNSS et de réanalyse de la station principale avec des données similaires provenant d'une ou plusieurs stations voisines. Chaque paire "station principale-station voisine" comprend quatre séries de base à partir desquelles six séries de différences sont formées. On propose et développe alors la stratégie suivante: chaque rupture détectée dans la série principale, un test de significativité d'un saut de moyenne au même instant est effectué dans les autres séries de différences formées, puis une règle statistique de prédiction permet d'attribuer la rupture testée au GNSS ou à la réanalyse. Lorsque plusieurs stations proches sont disponibles, une solution combinée est proposée.

Les nouveautés de notre approche résident dans (1) la prise en compte d'une dépendance temporelle évi-
dente et bien connue pour le test de significativité du saut à l'aide d'une méthode d'estimation GLS dans un
modèle de régression, et (2) du développement d'une règle de prédiction basé sur un algorithme supervisé connu.

Pour (1), la méthode consiste à caractériser dans un premier temps l'hétéroscédasticité et l'autocorrélation présentes dans le bruit des six séries de différences. Les modèles de bruit testés sont le bruit blanc, AR(1), MA(1) et ARMA(1,1). Il est constaté que seul un petit ensemble de séries se conforme au modèle de bruit blanc. Par conséquent, prendre en compte cette dépendance, une approche GLS, connue pour fournir des estimateurs statistiquement plus efficaces, est employée afin d'évaluer la signification des sauts. Le test se fait dans un modèle de régression prenant en compte l'hétéroscédasticité et le modèle de dépendance alors identifiée du bruit mais aussi le biais périodique dans la moyenne.

Pour (2), une règle prédictive est ensuite construite en formant un classificateur pour prédire l'origine d'une rupture à partir des résultats des tests obtenus sur un ensemble de jeu de données réelles. Deux difficultés principales se sont présentées: la réponse n'est pas connue sur les données réelles et le jeu de données est de faible taille et n'est clairement pas exhaustif pour espérer construire une bonne règle de prédiction. Pour relever ces défis, nous proposons de générer un ensemble de données synthétiques basé sur les résultats des tests des données réelles en utilisant une technique de rééchantillonnage.

Lorsque la méthode est appliquée aux 81 stations GNSS du jeu de données CODE REPRO2015, avec la
réanalyse ERA5 utilisée comme référence et le jeu de données NGL GNSS pour les stations voisines, 62%
des ruptures testées sont attribuées au GNSS, 19% à la réanalyse et 10% à des changements à la fois dans le
GNSS et la réanalyse. La majorité des points de changement identifiés est attribuée au GNSS ; cela constitue
un résultat favorable et est conforme aux attentes. Les résultats de la prédiction sont sensibles à la stratégie de
rééchantillonnage. Nous avons testé plusieurs variantes, en utilisant différents niveaux de signification pour les
tests, des échantillons équilibrés ou déséquilibrés et différentes règles d'agrégation lorsque plusieurs stations
proches sont disponibles. Les résultats des tests varient également en fonction du niveau de bruit, qui est
corrélé à la distance spatiale entre la station principale et les stations voisines. Les niveaux de bruit augmentent
considérablement avec la distance entre la station principale et les stations voisines dans un réseau peu dense. Ce
bruit élevé peut entraver les tests de signification dans les quatre séries non colocalisées, ce qui a inévitablement
un impact sur le résultat final de la prédiction.

## Conclusions et perspectives

En résumé, cette thèse contribue à la résolution de deux problèmes. Premièrement, elle a permis de mieux comprendre la sensibilité de la méthode de segmentation GNSSseg et des estimations de tendances aux propriétés des données GNSS et de référence (réanalyse). Deuxièmement, une méthode innovante d'attribution automatique a été développée, ce qui constitue une étape cruciale, avec la segmentation, au sein du processus global d'homogénéisation.

Plusieurs éléments de notre approche d'homogénéisation pourraient être améliorés et présenter des opportunités pour de futures recherches. Il s'agit notamment des éléments suivants:

1. L'approche d'attribution pourrait bénéficier d'un ensemble de données GNSS plus conséquent et plus exhaustif provenant de NGL repro3. Il serait possible d'affiner les critères de sélection des stations voisines, notamment en privilégiant des séries plus proches et avec moins de lacunes. Cela réduira sans aucune doute l'ampleur du bruit dans les séries de différences et améliorera la puissance du test du saut de moyenne. En outre, la taille plus importante de l'ensemble de données pourra permettre d'établir une nouvelle règle de classification plus performante, car si le nombre de changements testés est plus élevé, il y aura moins de répétitions dans l'ensemble de données synthétique construit par rééchantillonnage.

2. Un package R de la méthode d'attribution sera développé, permettant son accessibilité à l'ensemble de la communauté scientifique.

3. Avec le jeu de données NGL repro3, il sera aussi possible de tester d'autres approches de la méthode de segmentation mise en œuvre actuellement, comme l'utilisation de stations GNSS voisines comme référence au lieu d'une réanalyse. Lorsque la station principale et les stations voisines sont traitées de manière cohérente, les inhomogénéités restantes proviendraient principalement de changements spécifiques à la station, tels que des changements instrumentaux ou des changements dans l'environnement. Par conséquent, cette approche peut faciliter l'attribution des ruptures et conduire à des améliorations significatives dans les séries homogénéisées.

Dans le contexte du traitement des données GNSS, la segmentation s'avère très sensible aux petites variations dans la moyenne de la série. Cette sensibilité est influencée par plusieurs aspects clés du traitement des données GNSS, avec un accent particulier sur des facteurs tels que: la correction ZHD a priori, le modèle d'étalonnage de l'antenne/radome, la fonction de cartographie et l'angle de coupure de l'élévation. La cohérence et la fiabilité des produits GNSS générés par différents logiciels et incorporant diverses caractéristiques de traitement peuvent faire l'objet d'un examen minutieux.

Dans une perspective de recherche à long terme, la réalisation d'une analyse des tendances homogénéisées de la valeur de référence du GNSS offre une opportunité convaincante, en particulier dans le contexte des études sur le changement climatique. Cette analyse apporte des preuves empiriques pour répondre à une question clé: Comment la vapeur d'eau évoluera-t-elle dans un climat plus chaud ? La rétroaction de la vapeur d'eau est souvent estimée à l'aide d'un modèle climatique global (MCG), qui dépend fortement de la paramétrisation. Les jeux de données GNSS homogénéisées peuvent être comparés et utilisés comme référence pour le réglage des
paramètres dans les modèles climatiques, ce qui permet d'obtenir les "meilleures" estimations de la rétroaction de la vapeur d'eau et de ses incertitudes.
# Glossary

**ACF** Autocorrelation function

**AR** Autoregressive model

**ARMA** Autoregressive Moving Average model

**CART** Classification and Regression Tree

**CMIP** Coupled Model Intercomparison Project

**CODE** Center for Orbit Determination in Europe

**ERA-Interim** ECMWF Reanalysis - Interim

**ERA5** ECMWF Reanalysis version 5

**GNSS** Global Navigation Satellite System

**IERS** International Earth Rotation and Reference Systems Service

**IGS** International GNSS Service

**IPCC** Intergovernmental Panel on Climate Change

**IWV** Integrated Water Vapor

**JPL/NASA** NASA's Jet Propulsion Laboratory

**k-NN** k Nearest Neighbors

**LDA** Linear Discriminant Analysis

**MA** Moving Average model

**MERRA-2** Modern-Era Retrospective analysis for Research and Applications Version 2

**NGL** Nevada Geodetic Laboratory

**PACF** Partial autocorrelation function

**PW** Precipitable Water

**RF** Random Forest

**ZHD** Zenith Hydrostatic Delay

**ZTD** Zenith Tropospheric Delay

**ZWD** Zenith Wet Delay
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>1.1</td><td>Background and context</td><td>1</td></tr><tr><td>1.1.1</td><td>The role of water vapour in a changing climate</td><td>1</td></tr><tr><td>1.1.2</td><td>Consistency and homogeneity of observations, reanalyses and climate models</td><td>5</td></tr><tr><td>1.1.3</td><td>Inhomogeneities in GNSS IWV time series</td><td>7</td></tr><tr><td>1.2</td><td>Review of existing homogenization methods</td><td>9</td></tr><tr><td>1.3</td><td>Objectives and outline of this thesis</td><td>11</td></tr><tr><td>2</td><td>Data and mathematical tools</td><td>13</td></tr><tr><td>2.1</td><td>Data</td><td>13</td></tr><tr><td>2.1.1</td><td>GNSS data</td><td>13</td></tr><tr><td>2.1.1.1</td><td>GNSS observations and data processing methods</td><td>13</td></tr><tr><td>2.1.1.2</td><td>Modelling the tropospheric delays</td><td>15</td></tr><tr><td>2.1.1.3</td><td>Conversion of tropospheric delay to IWV</td><td>17</td></tr><tr><td>2.1.1.4</td><td>Important data processing aspects</td><td>18</td></tr><tr><td>2.1.1.5</td><td>GNSS tropospheric delay products used in this work</td><td>19</td></tr><tr><td>2.1.2</td><td>Reanalysis data</td><td>21</td></tr><tr><td>2.2</td><td>Mathematical tools</td><td>22</td></tr><tr><td>2.2.1</td><td>GNSSseg segmentation method</td><td>22</td></tr><tr><td>2.2.1.1</td><td>Model</td><td>22</td></tr><tr><td>2.2.1.2</td><td>Estimation procedure</td><td>23</td></tr><tr><td>2.2.1.3</td><td>Important features</td><td>24</td></tr><tr><td>2.2.2</td><td>Linear Regression</td><td>25</td></tr><tr><td>2.2.3</td><td>Linear stochastic models</td><td>27</td></tr><tr><td>2.2.3.1</td><td>Definition</td><td>27</td></tr><tr><td>2.2.3.2</td><td>Important features</td><td>29</td></tr><tr><td>2.2.3.3</td><td>Fundamental models</td><td>29</td></tr><tr><td>2.2.3.4</td><td>Note on the PACF</td><td>32</td></tr><tr><td>2.2.3.5</td><td>Note on spectral properties of stationary models</td><td>34</td></tr><tr><td>2.2.3.6</td><td>Combination of models</td><td>35</td></tr><tr><td>2.2.3.7</td><td>Model identification</td><td>36</td></tr><tr><td>2.2.3.8</td><td>Parameter estimation</td><td>36</td></tr></table>
<table><tr><td>2.2.4</td><td>Classification</td><td>37</td></tr><tr><td>2.2.4.1</td><td>The popular learning algorithms</td><td>38</td></tr><tr><td>2.2.4.2</td><td>Resampling methods</td><td>43</td></tr><tr><td>2.2.4.3</td><td>Imbalance class problem</td><td>45</td></tr></table>

3 Sensitivity of Change-Point Detection and Trend Estimates to GNSS IWV Time Series Properties 46

<table><tr><td>3.1</td><td>Summary</td><td>46</td></tr><tr><td>3.2</td><td>Paper No. 1</td><td>47</td></tr><tr><td>3.2.1</td><td>Abstract</td><td>48</td></tr><tr><td>3.2.2</td><td>Introduction</td><td>48</td></tr><tr><td>3.2.3</td><td>Materials and Methods</td><td>50</td></tr><tr><td></td><td>GNSS IWV Data</td><td>50</td></tr><tr><td></td><td>Reference IWV Data</td><td>54</td></tr><tr><td></td><td>Homogenization Method</td><td>55</td></tr><tr><td></td><td>Trend Estimation Method</td><td>56</td></tr><tr><td>3.2.4</td><td>Results</td><td>57</td></tr><tr><td></td><td>Segmentation Results</td><td>57</td></tr><tr><td></td><td>IWV Trend Estimates</td><td>70</td></tr><tr><td>3.2.5</td><td>Conclusions</td><td>75</td></tr></table>

4 Development of the attribution method 80

<table><tr><td>4.1</td><td>Paper No. 2: A statistical method for the attribution of change-points in segmented IWV difference time series</td><td>80</td></tr><tr><td>4.1.1</td><td>Abstract</td><td>80</td></tr><tr><td>4.1.2</td><td>Introduction</td><td>81</td></tr><tr><td>4.1.3</td><td>Data characterization</td><td>85</td></tr><tr><td>4.1.3.1</td><td>Data sets</td><td>85</td></tr><tr><td>4.1.3.2</td><td>Pre-processing</td><td>86</td></tr><tr><td>4.1.3.3</td><td>Characterization of the data and model building</td><td>86</td></tr><tr><td>4.1.4</td><td>Proposed tests for a fixed change-point</td><td>90</td></tr><tr><td>4.1.4.1</td><td>Regression model and different tests</td><td>90</td></tr><tr><td>4.1.4.2</td><td>Evaluation based on simulations</td><td>92</td></tr><tr><td>4.1.4.3</td><td>Application to real data</td><td>94</td></tr><tr><td>4.1.5</td><td>Predictive rule</td><td>95</td></tr><tr><td>4.1.5.1</td><td>Preliminary considerations</td><td>96</td></tr><tr><td>4.1.5.2</td><td>The proposed Cross-Validation Bootstrap (CVB) procedure</td><td>97</td></tr><tr><td>4.1.5.3</td><td>Application to the real data set</td><td>98</td></tr><tr><td>4.1.6</td><td>Conclusions and perspectives</td><td>101</td></tr><tr><td>4.1.7</td><td>Appendix</td><td>102</td></tr><tr><td>4.1.7.1</td><td>Test table</td><td>102</td></tr><tr><td>4.1.7.2</td><td>Computation of prior probabilities</td><td>103</td></tr><tr><td>4.2</td><td>Additional studies</td><td>104</td></tr></table>
<table><tbody><tr><td>4.2.1</td><td>Assessment of the stochastic model identification and parameter estimation</td><td>104</td></tr><tr><td>4.2.1.1</td><td>Simulation set-up</td><td>104</td></tr><tr><td>4.2.1.2</td><td>Noise model identification results</td><td>104</td></tr><tr><td>4.2.1.3</td><td>Parameter estimation results</td><td>112</td></tr><tr><td>4.2.1.4</td><td>Impact of gaps in the time series</td><td>113</td></tr><tr><td>4.2.1.5</td><td>Interpretation of results from real data</td><td>114</td></tr><tr><td>4.2.2</td><td>Understanding of test and prediction results</td><td>116</td></tr><tr><td>4.2.3</td><td>Strategies for further improving the prediction results</td><td>126</td></tr><tr><td>4.2.3.1</td><td>Correction of distance bias</td><td>126</td></tr><tr><td>4.2.3.2</td><td>Multiple testing correction</td><td>127</td></tr><tr><td>4.2.4</td><td>Conclusions</td><td>128</td></tr><tr><td><b>5</b></td><td><b>Conclusions and perspectives</b></td><td><b>130</b></td></tr><tr><td>5.1</td><td>Conclusions</td><td>130</td></tr><tr><td>5.2</td><td>Perspectives</td><td>132</td></tr><tr><td><b>A</b></td><td><b>Appendix: Preliminary results with an extended dataset</b></td><td><b>135</b></td></tr></tbody></table>
# Chapter 1

# Introduction

## 1.1 Background and context

### 1.1.1 The role of water vapour in a changing climate

Water vapour plays a central role in Earth's energy and water cycle in a variety of ways. Being relatively transparent to the (shortwave) solar radiation and mainly opaque to the (longwave) terrestrial radiation, it is the strongest greenhouse gas in the atmosphere and contributes to nearly 75% of the total greenhouse effect (Bengtsson, 2010). It is also a key actor of energy transport, by absorbing latent heat during evaporation over the oceans and releasing it through condensation in clouds. This exchange of latent heat is pivotal in driving atmospheric circulation at various scales, from moist convection in individual thunderstorm cells to the long-range planetary circulation and moisture transport (Schneider et al., 2010; Sherwood et al., 2010). The strength of moist processes such as convection and precipitation is tightly controlled by the amount of water vapour in the atmosphere, where high amounts of moisture can lead to extreme events (e.g. heavy precipitation and floods). Another fundamental property of water vapour is its short residence time in the atmosphere of about 8 days (Trenberth, 1998), which makes it extremely variable, both spatially and temporally, and thus challenging to observe and simulate in atmospheric models.

One of the fundamental controls of water vapour in the climate system is the Clausius-Clapeyron (CC) law which relates saturation vapour pressure $e_s$ to temperature, $T$ (Peixoto and Oort, 1996):

$$ \frac{de_s}{dT} = \frac{Le_s}{R_v T^2}. \qquad (1.1) $$

where $R_v$ is the gas constant for water vapour (461 J K$^{-1}$ kg$^{-1}$), and $L$ represents the latent heat of vapourisation ($2.5 \times 10^6$ J kg$^{-1}$). If we integrate (1.1) between $T_0$ and $T$, we obtain:

$$ e_s(T) = e_0 \exp\left\{-\frac{L}{R_v}\left(\frac{1}{T} - \frac{1}{T_0}\right)\right\}. \qquad (1.2) $$

where $e_0 = e_s(T_0)$.

Because temperature decreases at a mean rate of ~ 6.5 K km$^{-1}$ in the troposphere, the saturation vapour pressure varies roughly by 4 orders of magnitude between the surface and the tropopause. This steep vertical
gradient plays a major role in convection. Moisture in a raising air parcel, even if not initially saturated, can thus quickly reach the saturation level above which the latent heat released by condensation will enhance its buoyancy and may ultimately lead to deep convection. This mechanism is a major humidification process of the upper troposphere and the lower stratosphere.

The actual amount of water vapour in the atmosphere can be quantified by different variables. Relative humidity (RH) is the quantity most used for general applications (Peixoto and Oort, 1996). It is defined as $U = e/e_s$, the ratio of vapour pressure, $e$, of an individual air parcel at a temperature $T$, to saturation vapour pressure of air at the same temperature, $e_s(T)$. A collection of saturation vapour pressure formula, $e_s$ as a function of $T$, can be found here (http://cires1.colorado.edu/voemel/vp.html). Specific humidity is another widely used quantity in meteorology which expresses the fraction of mass of water vapour in an air parcel at temperature $T$ and pressure $P$ to its total mass: $q = \rho_v/\rho$, where $\rho_v = e/R_vT$ and $\rho = P/RT$ are the volumic masses (densities) of water vapour and the mixture of dry and air and water vapour, respectively. Relative humidity differs both qualitatively and quantitatively from other moisture variables (Peixoto and Oort, 1996). The distribution of humidity and temperature as a function of height are fundamental for many atmospheric phenomena, e.g. initiation of convection. They can be measured in-situ by balloon-borne sensors (radiosondes). Because it is constrained by the Clausius-Clapeyron law, specific humidity is quickly decreasing with height, i.e. most of the water vapour in an atmospheric column is concentrated in the lower troposphere (about 90% lies below 5 km). For this reason, the Total Column Water Vapour (TCWV), also referred to as Integrated Water Vapor (IWV) or Precipitable Water (PW), is another relevant variable often used to quantify the amount of water vapour in the atmosphere. TCWV and IWV are defined as the vertical integral of water vapour density in the atmospheric column:

$$ \text{IWV} = \int_{0}^{\infty} \rho_{v}(z) dz. \qquad (1.3) $$

where $z = 0$ refers to the surface and $z = \infty$ to the top of atmosphere. IWV is expressed in units of kg m⁻² and represents the total mass of water vapour in an unit-area column of atmosphere. A typical mi-latitude summer value would be IWV = 25 kg m⁻². PW is defined as the integral of $\rho_v(z)/\rho_l$, where $\rho_l=1000 \text{ kg m}^{-3}$ is the density of liquid water. It is usually expressed in mm or cm and can be interpreted as the equivalent height of liquid water if all the water vapour in an unit-area column of atmosphere would be condensed and precipitated. IWV is a quantity that can be measured by Global Navigation Satellite System (GNSS) and microwave radiometers.

Figure 1.1, borrowed from Parracho et al. (2018), shows the geographical distribution of mean IWV in boreal winter and summer from ERA-Interim reanalysis and GNSS stations. IWV reaches its maximum around the equator, ranging from about 20 to 60 kg m⁻², due to intense evaporation over the oceans caused by warm sea surface temperatures and strong surface winds. At higher latitudes and over the continents, IWV decreases because of limited moisture supply over land masses and lower water holding capacity of the cooler atmosphere. Another noteworthy observation is that the reanalysis is able to captures strong spatial gradients in good agreement with GNSS data.

In the context of climate change, water vapour plays a significant role in enhancing global warming through a positive feedback mechanism (Bony et al., 2006; Colman and Soden, 2021). This feedback mechanism acts in
<figure><img src="image_2.png" /></figure>

Figure 1.1 – Mean value of IWV from ERA-Interim between 1995 and 2010 for December-January-February (DJF) (a) and June-July-August (JJA) (b). Filled circles correspond to IWV retrieved by Global Positioning System (Parracho et al., 2018).

response to the forcing effect of other greenhouse gases, such as CO2. Climate model simulations have shown that, in the scenario of doubling CO2, temperature increase is amplified by a factor of 2 to 3 due to the water vapour feedback compared to a "no-feedback" scenario (Manabe and Wetherald, 1967; Held and Soden, 2000). The mechanism is the following: an increase in non-condensable, long-lived, greenhouse gases (e.g. CO2, CH4 or NO2), induces a rise in temperature through greenhouse radiative forcing. Guided by the Clausius-Clapeyron (CC) law, this temperature increase results in a proportional elevation of saturation vapour pressure. Given that water vapour itself is a greenhouse gas, an increase in atmospheric water vapour enhances the absorption of both longwave and shortwave radiation while reducing outgoing terrestrial radiation. This, in turn, amplifies the warming effect. Moreover, the intensification of the water cycle, including heavy precipitation and flood events, is closely linked to the increase in low-altitude water vapour (Douville et al., 2021).

A key question is, what actual changes in humidity and temperature can we expect in a warming climate? (Colman and Soden, 2021) Global mean surface temperature in the first two decades of the 21st century (2001–2020) was 0.99 [0.84 to 1.10] °C higher than 1850–1900 (Gulev et al., 2021, p. 320). Guided by the CC law, the corresponding rate of increase in saturation vapour pressure at T=273 K is ~ 7% K⁻¹, but how much is the actual water vapour changing?

Both observations and climate models have consistently supported the assumption that relative humidity remains relatively close-to-unchanged on a global scale in response to warming (Held and Soden, 2006; Trenberth et al.,
2013; Hartmann et al., 2013; Colman and Soden, 2021; Douville et al., 2022). Unchanged RH implies an increase in specific humidity and in IWV, which was also confirmed in IWV observations from various data sources (Hartmann et al., 2013; Bock et al., 2014; Wang et al., 2016; Parracho et al., 2018).

Figures 1.2 and 1.3 show the global mean (90°S-90°N) lower-tropospheric temperature (LTT) anomalies and the near global (60°S-60°N) IWV anomalies over the past decades from observations and reanalyses (Dunn et al., 2023). Both temperature and IWV show a long-term tendency towards higher values, superposed to quite large interannual to decadal variability. The linear trend estimates of the various temperature data sets of Figure 1.2 are in the range 0.17-0.22 K decade⁻¹ for the period 1979-2022 (Dunn et al., 2023, p. S36). The trends from the three reanalyses (ERA5, MERRA-2, and JRA55) are in the range 0.13-0.23 K decade⁻¹ for the more recent period, 1995-2021. The corresponding IWV trends are ranging from 1.3 to 1.6% decade⁻¹ and the CC scaling factors take values from 6.9 to 7.9% K⁻¹ [O. Bock, personal communication]. They align closely to the 7% K⁻¹ predicted by the Clausius-Clapeyron equation at T = 273 K, and thus confirm that the global mean water vapour has increased at nearly constant RH over the past 2.5 decades.

<figure><img src="image_3.png" /><figcaption>Figure 1.2 – Monthly average global lower-tropospheric temperature (LTT) anomalies (°C; 1991–2020 base period) for (a) radiosonde, (b) satellite, and (c) reanalysis datasets. Time series are smoothed using a 12-month running average. Annual averages are displayed for the Radiosonde Atmospheric Temperature Products for Assessing Climate (RATPAC) dataset (Fig. A2.7 from Dunn et al., 2023).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 1.3 – Global mean total column water vapour annual anomalies (kg m⁻²) over (a) land and ocean, (b) ocean only, and (c) land only from observations and reanalyses (ERA5, MERRA-2, JRA-55). The shorter time series from the observations have been adjusted so that there is zero mean difference relative to the ERA5 results during their respective periods of record (Fig. 2.25 from Dunn et al., 2023).</figcaption></figure>
### 1.1.2 Consistency and homogeneity of observations, reanalyses and climate models

While a global upward trend in IWV has been deemed very likely, it is essential to evaluate the level of uncertainty in the trend estimates and understand the origin of discrepancies between the various data sources. Figure 1.3 reveals differences between the various observational techniques, as well as between the various reanalyses. The latter depart from each other especially in years before 1995 where the number of assimilated data is much lower than in more recent years. Changes in the geographical distribution and instrumentation of stations over time are a major source of spatial and temporal inhomogeneity for reanalyses (Hersbach et al., 2020) and their capacity for estimating decadal trends has been questioned (Thorne and Vose, 2010). Nevertheless, reanalysis data are extensively used to complement observational records for the evaluation of climate model simulations (Boucher et al., 2020).

In the framework of the IPCC Fifth Assessment Report (AR5), a large dispersion in IWV and lower tropospheric temperature (LTT) trends was evidenced among the Fifth Coupled Model Intercomparison Project (CMIP5) climate model simulations in tropical regions (20°S-20°N). Although all simulated trends align well on a 5.9% K⁻¹ CC line (fairly consistent with CC and the assumption of constant RH), the magnitudes of their LTT and IWV trends ranged from 0.1 to 0.5 K decade⁻¹ and 0.5 to 2.5% decade⁻¹, respectively (Flato et al., 2013, Fig. 9.9, p. 774). On the opposite, the trends associated with MERRA-2 and ERA-Interim reanalyses, and with two satellite products, are relatively far from the CC line, and thus rather suggest a long-term change in relative humidity. It is not known whether these discrepancies are due to remaining inhomogeneity in the observational data and/or reanalysis results, or due to problems with the climate simulations. Similar discrepancies have been reported with the more recent CMIP6 model simulations, ERA5 reanalysis, and satellite observations by Santer et al. (2021).

In a study by Allan et al. (2022), a comparison was drawn among satellite observations, ERA5 reanalysis, and CMIP6 climate model simulations. For the model simulations, they distinguished the *historical* (coupled) simulations and *amip* (atmosphere-only simulation) experiments for a subset of models. They found median IWV trend estimates of 1.9 and 1.1% decade⁻¹ for the two types of simulations (period 1988-2014) compared to 0.78% decade⁻¹ for ERA5 and 1.0% decade⁻¹ for observations which combined infrared and microwave satellite data and ground-based humidity observations from stations. The *amip* simulations exhibited quite consistent trends with the observations. In contrast, CC scaling factors of the *historical* simulations (5.6% K⁻¹) showed better agreement with observations (5.5% K⁻¹) and ERA5 (5.8% K⁻¹) than the *amip* simulations (4.8% K⁻¹).

More recently, Douville et al. (2022) investigated the feasibility of constraining climate model projections of global mean IWV with the help of observations and reanalysis data using a Bayesian statistical method. Constraining the forced response of historical CMIP model simulations with global mean surface temperature observations since 1850 and global mean IWV data from GNSS observations and reanalyses after 1994, they could reduce the spread in global mean IWV projected to the end of the 21st century by 39%. The 5-95% confidence interval of the constrained CC scaling factors was also reduced (6.5 to 7.6% K⁻¹), with a median close to 7% K⁻¹. The projections also indicate a strong geographical disparity in the water vapour change, with a multi-model mean exceeding 12 kg m⁻² in the tropics or 50% in desertic and polar regions, which may be
connected to deviations from CC-scaling, especially over land where constant RH cannot be maintained due to lack moisture supply from the surface despite the increasing evaporative demand. In a similar study using surface RH measurements, constrained model projections show an inevitable continental drying, especially in the northern midlatitudes (Douville and Willett, 2023).

<figure><img src="image_5.png" /><figcaption>Figure 1.4 – Absolute (top) and relative (bottom) trend of the IWV. The map shows results from ERA5 (left) and MERRA-2 (right) reanalyses, while filled circles correspond to GNSS stations (courtesy O. Bock, 2023).</figcaption></figure>

Figure 1.4 presents both absolute and relative trends of monthly IWV anomalies obtained at 240 GNSS stations (Bock, 2020b), alongside the MERRA-2 and ERA5 reanalyses, using the same data set as for Figure 1.3. It is important to note that, although there is a general positive trend observed across all datasets, substantial variations in the magnitude of the trends is observed between the reanalyses in various regions, with sometimes opposite signs, especially over the maritime continent, but also in the Atlantic Ocean, South America, South Asia, and Antarctica. One may note that the trend expressed in kg m⁻² decade⁻¹ is very close to zero in Antarctica because the mean IWV is usually very low there. However, the relative trend (in % decade⁻¹) shows both positive and negative values, but with relatively poor consensus between ERA5 and MERRA-2. Parracho et al. (2018) also noted pronounced disparities between ERA-Interim and MERRA-2 over Antarctica, except at the location of the GNSS stations. Although GNSS data are not assimilated in these reanalyses, these GNSS stations are located in places where other meteorological measurements are collected and are probably assimilated.

Discrepancies between GNSS and the ERA-Interim at a number of stations were reported by Parracho et al. (2018). Inspection of the IWV time series for these stations revealed inhomogeneities in the GNSS data, in
the form of jumps and drifts, which were identified as instrumentation malfunctioning and/or replacement. Compared to Parracho et al. (2018)'s results, Figure 1.4 shows much better agreement between reanalyses and GNSS (where reanalyses agree), suggesting that the updated GNSS product used in Figure here is more homogeneous than the older one (this point will be confirmed in Chapter 3. However, inhomogeneities in the reanalyses cannot be excluded. For example, (Schroeder et al., 2016) detected jumps in IWV from reanalyses and satellite observations, and noted that several change-points in the reanalyses coincided with changes in the observing system (e.g. start and end of assimilation of satellite data in the reanalyses).

### 1.1.3 Inhomogeneities in GNSS IWV time series

In the context of atmospheric water vapour monitoring, GNSS offers reliable observational data with a high temporal resolution available under all weather conditions (Guerova et al., 2016). As long as these observations are not assimilated in reanalyses, they remain useful independent data for evaluating the quality of the latter, as well as climate model outputs (Parracho, 2017). Nevertheless, the homogeneity of GNSS data is susceptible to various factors, including instrumental changes, processing variations, and environmental influences (Vey et al., 2009; Bock et al., 2010; Ning et al., 2016). Jumps in IWV series of 1 to 2 kg m⁻² have been observed which are susceptible to obscure most of the underlying climatic trends which are in the range of -2 to +2 kg m⁻² decade⁻¹ at stations (Figure 1.4). Homogenization of GNSS IWV series appeared as a necessity for climate applications and the first consistent effort to this end was undertaken during the COST GNSS4SWEC project (Bock O. Pacione R. Ahmed F. Araszkiewicz A. BaŁdysz et al., 2020). Several available tools used by the climate community for the homogenization of temperature and precipitation data were tested in a benchmarking exercise organized in the framework of this project (Van Malderen et al., 2020). At the same time, Quarello (2020) developed a specialized method for the detection of change-points in GNSS minus reanalysis IWV difference series, which proved to be one of the best tested by Van Malderen et al. (2020). This thesis builds on the detection tool proposed by Quarello (2020) to include it in a full homogenization method. Before going into further details of this method, we explain the reason why GNSS minus reanalysis difference series need to be considered.

The top plot in Figure 1.5 illustrates the daily IWV GNSS series recorded at the ALIC station in Alice Springs, Australia. The IWV series at this station displays a pronounced temporal variation, with values ranging from under 5 kg m⁻² in austral winter to over 50 kg m⁻² in summer. Such substantial seasonal variability can pose challenges in detecting small jumps in the raw IWV time series. The same is observed with temperature time series and other climatic variables (Easterling and Peterson, 1995; Mitchell and Jones, 2005). To address this issue, following the "relative segmentation" approach developed by climatologists, we differentiate the GNSS IWV series with respect to a reference series. The lower plot in Figure 1.5 thus shows the IWV difference series between the GNSS and ERA5 reanalysis data extracted at the location of the station. The strong seasonal variation seen in the raw IWV series is clearly reduced, thereby improving the detectability of small inhomogeneities. The figure also presents the segmentation results (change-points are represented by vertical dashed lines) obtained with the GNSSseg method developed by Quarello (2020) (the mathematical aspects of this method will be explained in subsection 2.2.1). However, it is important to notice that, although significantly reduced, the seasonal variation in the difference series is not entirely removed. Additionally, the day-to-day scatter the difference series also exhibits strong seasonality. Both features are due to representativeness differences between the GNSS point observations and the model grid cells (Bock and Parracho, 2019), and are susceptible to induce false change-
point detections (Bock et al., 2019). In Quarello (2020)'s segmentation method, these features are included in the mathematical model in the form of a periodic bias and a monthly variance, respectively. The estimates of these two features are shown in Figure 1.5 as the purple line for the bias and the cyan line for the monthly variance.

At station ALIC, the segmentation tool detected five change-points, marked by the vertical dashed lines. They delimit changes in the mean, drawn as the broken red line, which are in the range 0.5 - 1.5 kg m⁻². Since the segmentation is conducted on the difference series, these inhomogeneities could originate from either the GNSS or the ERA5 data. Metadata documenting instrumental changes at the GNSS station can serve as a valuable source for validating the detected change-points. This information is reported on Figure 1.5 for the ALIC stations. Small colored triangles indicate the time position of instrumental (receiver, antenna, and/or radome) changes as well as the position of a change in the processing procedure (blue triangle). Among the five detected change-points, two coincide well with equipment changes: a receiver change in 2000 and a receiver plus antenna change in 2011. However, the other change-points are a bit far from known equipment changes to be surely attributed to GNSS origin.

<figure><img src="image_6.png" /></figure>

Figure 1.5 – Daily time series of IWV GNSS (top) and IWV difference between GNSS and ERA5 reanalysis (bottom) at station ALIC (Alice Springs, Australia). The bottom plot includes segmentation results superposed (the red line shows the mean for each segment, the vertical black dashed lines indicate the change-point positions, the cyan and purple lines at the bottom represent the estimated monthly mean variance and the functional modeling the periodic bias, respectively). Colored triangles on top of each plot indicate the position of known instrumental or processing changes from the GNSS metadata.
### 1.2 Review of existing homogenization methods

Homogenization refers to the process of identifying and correcting (artificial) inhomogeneities within data series. In this study, our focus is on the detection of abrupt changes or change-points within the data, as they can significantly impact the trend estimates. Two approaches have been commonly used for the detection step: either directly on the raw series, known as the "absolute" approach, or on the series of differences (as mentioned in the previous subsection), known as the "relative" approach. To mitigate the impact of large seasonal variations in the raw signal, the absolute approach needs either to include a mathematical model for the seasonal variations or to work on anomalies (differences with respect to the mean annual variation) instead of the raw signal. In the relative approach, this problem is greatly reduced but the introduction of a reference series also introduces uncertainty on the origin of the detected change-points if the reference series is not perfectly homogeneous. Both approaches have been discussed and an exhaustive number of existing software packages have been assessed in the COST Action HOME (Venema et al., 2012). Their conclusion was unequivocal on the higher performance of the relative approach, although the attribution step remains a limiting factor.

In the rest of this section, we describe the different steps necessary to the relative homogenization approach.

**Data preparation.** The first task is to compute the series of differences between the target series and a reference series. The reference series can be taken from a single nearby station or computed as a composite from several nearby stations. It can also be taken from a different data source, such as reanalysis, as long as it shares a common climate signal with the target series (the goal is to suppress this signal by taking the difference of series). This step may actually consist of several sub-steps to minimize the contribution of various error sources present in the raw series or introduced in the series of differences, such as outliers and representativeness differences. The representativeness differences may usually recover two types of problems. The first one was already mentioned above for the case of the reanalysis. Indeed, reanalysis data are representative of aerial averages on the order of (twice) the product grid spacing (e.g. about 50 km for the 0.25° grid size of ERA5), while GNSS observations are more representative of the atmosphere in the close vicinity of the antenna, or a few km around it. The second type arises when a nearby station is used as a reference and there is a significant climatic difference between stations (e.g. when the distance is too far, e.g. hundreds of km, or even a shorter distances of a few km in mountainous areas). Another important feature is to make the IWV data consistent in the vertical dimension because water vapour density and thus IWV are quickly decreasing with altitude. Therefore, a correction of the IWV data due to vertical distance is required (Bock et al., 2022).

**Segmentation.** The central step of homogenization the segmentation which aim is to detect change-points (usually at unknown time). This can be done manually, by the visual inspection of time series, fully automatically, using a statistical method, or in a hybrid way where a statistical method is used to detect potential change-points which are validated afterwards by an expert. A comprehensive review of segmentation methods used for the homogeneization of climate data has been done in Quarello (2020). The list of the current segmentation methods, along with associated software (if available), can be found in Table 2.1 of that thesis, as well as in references therein. Most of these methods have been assessed in the COST Action HOME (Venema et al., 2012), for monthly temperature and precipitation observations, and more recently in COST Action GNSS4WSWEC for GNSS IWV observations (Van Malderen et al., 2020). These segmentation methods are classified based on their model type (parametric or non-parametric), approach (frequentist or Bayesian), infer-
ence procedure (likelihood, penalized likelihood, or test), and search algorithm (optimal or sub-optimal), as shown in Figure 2.7 in Quarello (2020). Most of these methods operate using the relative approach, which is a common method for reducing strong seasonal variations and increasing the detection power of the segmentation.

In this study, we employ the segmentation method developed by Quarello et al. (2022) in a parametric and frequentist framework. This method allows to detect all change-points simultaneously in an efficient manner or using an optimal search computationally through the use of dynamic programming. This method has been specifically developed for detecting change-points in Integrated Water Vapour (IWV) difference series between GNSS and reanalysis data. In particular, a change-point detection in the mean model is proposed taking into account the characteristics of the target series such as periodic bias and varying variance as shown in the example in Figure 1.5. The inference of the model is done through a two-step strategy: first, for a given number of change-points, the parameters of the model (the means of the segments, the variances of each month, the periodic function and the change-points) are estimated using the classical maximum likelihood procedure. Second, the number of change-points is estimated. This segmentation method effectively addresses two fundamental questions in the procedure, an algorithmical and a statistical issues:

★ Estimation of change-points position: the classical dynamic programming (DP) algorithm is used to locate change-points. This algorithm is now well known to retrieve the change-point efficiently (fastly and exactly) computationally speaking. This algorithm can be used if and only if the quantity to be optimized is additive with respect to the segment (see details in sub-section 2.2.1). The proposed inference method has been developed in order to be able to use it.

★ Estimation of the number of change-points: this issue is resolved using penalized criteria, as classically used in segmentation frameworks. A penalty term is added to the inference criterion (e.g., log-likelihood) to account for the model's complexity, i.e. the number parameter to be estimated.

Quarello (2020) extensively evaluated this segmentation method by comparing its performance using four penalized criteria, in both simulation and real data scenarios. The final version was released as the R package named "GNSSseg" and published (Quarello et al., 2022). Furthermore, this method underwent a comparison with other segmentation methods in a benchmark exercice coordinated by Van Malderen et al. (2014) where it demonstrates strong potential and revealed as one of the best methods.

**Attribution/validation.** Change-points identified through segmentation can originate from either the target or the reference series. Therefore, attributing the detections to the right series is crucial before correcting the raw IWV data for the jumps. There are two primary strategies to address this attribution problem. The first approach aims to improve the homogeneity of the reference series, thus making it possible to attribute all change-points to the target series. Typically, these reference series are generated by averaging data from multiple nearby stations (Alexandersson, 1986; Menne and Williams, 2005; Thorne and Vose, 2010). However, in practice, even the composited series often have non-negligible inhomogeneities. The second approach looks for the origin of the change-point through pairwise comparisons. This can be achieved semi-automatically, employing statistical inference and manual cross-checking with historical information (station metadata), in an iterative way (Caussinus and Mestre, 2004). However, there could be gaps or mistakes in the available metadata. Another method, proposed by Menne and Williams (2009), fully automates the process by assigning change-points to the station that has the most detections, when segmenting all possible combination series of
differences. This method therefore requires a large number of nearby stations to be effective. This approach
also leverages station metadata when available and attributes detected change-points to the closest known events
in the station's history within specified confidence limits.

**Correction.** The final step involves correcting the target time series for the jumps in mean corresponding to the attributed or validated change-points only. The traditional way of doing this is to subtract the estimated offsets from the GNSS IWV series while leaving the most recent homogeneous segment unchanged (the rationale for this is that the most recent data may have smaller absolute bias than the older data) (Van Malderen et al., 2020).

**1.3 Objectives and outline of this thesis**

The objective of this study is to establish a comprehensive homogenization procedure for the IWV time series
derived from GNSS data using the segmentation method developed by Quarello (2020). Therefore, this thesis
focuses on two major tasks:

★ The first task involves a rigorous investigation of the robustness, sensitivity, and inherent limitations of the segmentation method. This endeavor is particularly crucial for practical real-world applications. The scientific questioning and factors influencing this analysis encompass the quality of the GNSS data which primarily depend on the GNSS data processing methodologies (do more recent GNSS solutions have reduced inhomogenieties?), the length of the analysed time series (as the length is increasing with years, what is the impact on the segmentation performance?), variations in auxiliary data employed in IWV conversion, and the selection of the data source for the reference time series (multiple choices can be made among various reanalysis products). Additionally, the study explores the extent to which changes in the segmentation results impact the subsequent IWV trend estimates (is it more beneficial to detect and correct more small jumps or fewer but larger ones?).

★ The second task centers on the development of an automatic attribution method. The method is intended
to operate as a post-processing step separate from segmentation and should be applicable even in cases of
sparse GNSS networks where only a few nearby station are within a reasonable (100-200 km) distance.

The thesis is structured as follows: Chapter 2 provides an introduction to the datasets employed and to the
fundamental mathematical concepts and tools utilized in this study. This chapter covers: i) the retrieval of IWV
from GNSS observations, including the basic principles of GNSS data processing; ii) a brief description of the
mathematics of the GNSSseg segmentation method developed by Quarello (2020); iii) fundamentals of linear
regression methods and linear stochastic model identification, with a special focus on the mixed autoregressive
- moving average (ARMA) models that will be effectively utilized for modeling the noise in the IWV difference
series; and iv) an introduction to a number of classification algorithms employed in the development of the
attribution method. Chapter 3 is dedicated to the first task. It is presented as a paper published in the Atmosphere
journal. In this paper, it is important to note that the segmentation results and trend estimates did not yet include
the attribution step which was developed afterwards. Especially, for the trend estimates, the results include
two versions: one with all detected change-points corrected (as if only GNSS data were inhomogeneous) and
one with only the change-points that could be directly validated with the help of metadata. The main caveat
with the latter is that changes due to external factors, not reported in the metadata, are not taken into account.
Some true change-points may thus remain in the corrected time series. Chapter 4 comprises two main sections.
First is a paper, submitted to the International Journal of Climatology (presently under review) which described and applies the developed attribution method on a real data set of 81 main stations. Second in this chapter is additional research conducted alongside the paper's topic. Finally, Chapter 5 summarizes the main findings presented in this thesis and discusses future perspectives. The thesis contains also an Appendix section which presents some preliminary results of the application of the full homogenization process to an extended GNSS dataset comprising more than 6000 stations.
# Chapter 2

# Data and mathematical tools

## 2.1 Data

### 2.1.1 GNSS data

GNSS Integrated Water Vapor (IWV) derives from the propagation delay of radio waves in the atmosphere on their path from the GNSS satellites to the ground-based receivers (Figure 2.1). Data processing with scientific software and use of a number of auxiliary geodetic products is required to achieve estimates of propagation delays and subsequent IWV data with accuracy compatible for meteorology and climatology. In this sub-section we summarize the main steps and highlight the crucial aspects of data processing intervening in the elaboration of GNSS IWV data.

#### 2.1.1.1 GNSS observations and data processing methods

GNSS are satellite radio navigation systems providing global positioning everywhere on Earth (Hofmann-Wellenhof et al., 2007). Their principle is based on the transmission of radio signals by satellite constellations which allow receiver-equipped users to access precise three-dimensional positioning in real time with precision 1-10 m. The most popular GNSS is the US GPS (Global Positioning System), which is fully operational since 1994. It has been complemented over the years by the Russian GLONASS, the European Galileo, and the chinese Beidou, among others (Teunissen and Montenbruck, 2017a, part B). Ground segments, relying on tracking networks and data processing facilities, are associated to the space segment to provide essential navigation products such as satellite ephemerides, clock synchronization information, and Earth orientation parameters (the latter are necessary to relate the inertial reference frame, which is natural to satellite orbits, to the Earth-fixed frame used to express the receiver coordinates). Highly accurate versions of these products are computed in near real time by the International GNSS Service (IGS) and International Earth Rotation and Reference Systems Service (IERS) analysis centres for scientific usage.

Some of the inhomogeneities that we are tracking in the GNSS IWV time series can be related to changes in the observations, data processing procedure, and auxiliary products used for the processing and/or the conversion. To understand how these features can impact the propagation delays and IWV, it is necessary to provide some insight into the data processing procedure.
Geodetic positioning relies on the analysis of the so-called "phase observations". A single phase observation between a receiver *r* and satellite *s* can be formalized by the following equation (Bock, 2012, p. 6):

$$L_r^s = \rho_r^s + c(\Delta t_r - \Delta t^s) + \lambda N_r^s + \Delta \rho_{rel} - \Delta \rho_{iono} + \Delta \rho_{tropo} + \Delta \rho_{ant}^s + \Delta \rho_{ant,r} + \epsilon. \quad (2.1)$$

where:

* $\rho_r^s$ is the geometric distance between the receiver's antenna reference point (ARP) and satellite ARP,

* $c$ is the speed of light in a vacuum,

* $\Delta t_r$ and $\Delta t^s$ are the receiver and satellite clock offsets from the GNSS time reference,

* $\lambda$ is the wavelength of the radio signal,

* $N_r^s$ is the integer phase ambiguity,

* $\Delta \rho_{rel}$ account for relativistic effects,

* $\Delta \rho_{iono}$ is the propagation delay in the ionosphere,

* $\Delta \rho_{tropo}$ is the propagation delay in the troposphere,

* $\Delta \rho_{ant}^s$ and $\Delta \rho_{ant,r}$ are the phase delays in the satellite and the receiver antennas,

* $\epsilon$ accounts for other minor errors.

Basically, GNSS data processing is a regression problem in which unknown parameters of the observation equation are estimated in order to minimize the errors between the observed and predicted phase measurements. The IGS standards in "static" geodetic processing (i.e. positioning of ground-fixed antennas) consists in estimating the parameters in (2.1) in 24 hour batches or "sessions". Observations are typically sampled at a 30-sec rate, while the sampling of parameters is much lower (e.g. 1 set of receiver coordinates per session, 1 hour sampling for tropospheric delay parameters, etc.). However, not all parameters can be estimated simultaneously because of multi-colinearity. Namely, satellite positions and clock delays are usually fixed to their values computed by the IGS (so-called IGS satellite products), and relativistic effects which can be predicted by theory with a high degree of precision (Zhu and Groten, 1988) are corrected a priori. Receiver clock offsets can be either estimated "by epoch" (one parameter per observation) or eliminated by combining the simultaneous observations from two satellites and two receivers (so-called double differences) but they require that the data from multiple stations are used once. The ionospheric delays are usually eliminated by combining instantaneous observations from 2 or 3 frequencies, thanks to the dispersive nature of the refraction index of the ionospheric plasma region at the GNSS frequencies.

Phase delays in the satellite and receiver antennas are an important source of bias that needs be corrected (Hofmann-Wellenhof et al., 2007). Empirical models are used therefore, comprising a constant offset component and a variable hemispheric model which accounts for variations depending on the angle of incidence and frequency of the radio signal with respect to the antenna's orientation. These variations are contingent on the specific antenna model and type in use. For the satellite, the antenna offset is expressed relative to its center of mass and is on the order of 1 m. The offsets for the receiver antennas are given relative to the ARPs and are of order 0.1 meter. The variations are typically one order of magnitude smaller but need also to be corrected.
Eventually, station positions, receiver clock offsets, tropospheric delay, and ambiguity numbers are the main parameters that are estimated.

When the observations are close to zenith, several of these parameters are tightly correlated (station height, receiver clock offset, tropospheric delay, and ambiguity number) and biases in the a priori correction models (e.g. antenna phase model offsets) and/or in the auxiliary products can map onto all these parameters (Teunissen and Montenbruck, 2017b). One way to reduce this problem is to include observations at low elevation angles.

<figure><img src="image_7.png" /><figcaption>Figure 2.1 – Schematic presentation of individual slant delays (SD) from GNSS satellites and their mapping to zenith tropospheric delay (ZTD), (Guerova et al., 2016).</figcaption></figure>

Apart from the factors explicitly introduced in equation (2.1), the error term, $\epsilon$, includes all other error sources, namely the "multipath" effect, which is a spurious phase modulation induced by reflective surfaces in the vicinity of the receiving antenna. Changes in the reflective properties of the environment around the receiver's antenna can modify this effect and be a source of inhomogeneities in the estimated parameters (positions, tropospheric delays, etc.).

### 2.1.1.2 Modelling the tropospheric delays

GNSS signals undergo both retardation and bending as they cross the troposphere and the lower stratosphere, because of the varying index of refraction along their propagation path. The combined effect is referred to as the "total tropospheric delay". It is also sometimes called "excess path length", in reference to geometrical optics, and denoted $\Delta L$, but it is identical to the term $\Delta\rho_{tropo}$ introduced in equation (2.1). Below we use the traditional notation, following Bevis et al. (1992), which writes:

$$ \Delta L = \int_{L} n(s)ds - G, \qquad (2.2) $$

where $n(s)$ corresponds to the index of refraction of the air at the position $s$ along the curved trajectory of the ray path denoted as $L$. The parameter $G$ signifies the straight-line geometrical path length through the atmosphere
assuming $n = 1$. A modified form of Equation (2.2) writes:

$$ \Delta L = \int_{L} [n(s) - 1] ds + S - G. $$

where, $S = \int_{L} ds$, is the curved path length along $L$. This equation provides a breakdown of the excess path length, with the first term on the right side accounting for the retardation effect (reduction the wave velocity due to $n > 1$), $\Delta L_r$, and the second term representing the bending effect (a purely geometrical effect). The bending effect is much smaller than the retardation at most elevation angles above 0°. Negelecting the bending effect, the first term is often re-expressed in terms of atmospheric refractivity $N = 10^6 \times (n - 1)$ as follows (Bevis et al., 1992):

$$ \Delta L = 10^{-6} \int_{L} N(s) ds \qquad (2.3) $$

The refractivity of the air can be decomposed into two basic groups of constituents, dry air and water vapor molecules. Following (Davis et al., 1985), the tropospheric delay is reformulated to account for this distinction:

$$ \Delta L = 10^{-6} \int_{L} N_{h}(s) ds + 10^{-6} \int_{L} N_{w}(s) ds = \Delta L_{h} + \Delta L_{w}. \qquad (2.4) $$

where the first integral is referred to the "hydrostatic delay", primarily due to the delay induced by dry air, and the "wet delay", due to the water molecules. Furthermore, Saastamoinen (1973) first noted that the hydrostatic delay integral taken in the zenith direction is simply proportional to surface pressure (by use of the hydrostatic equation in the atmospheric column). This property offers a strong constraint on the determination of the total delay in Equation (2.1) as the Zenith Hydrostatic Delay (ZHD) can be corrected a priori from auxiliary surface pressure data (from measurements or atmospheric model outputs), leaving only the Zenith Wet Delay (ZWD) as unknown parameter. Because the phase observations are collected in practice from satellites in any direction (see Figure 2.1), mapping functions are used to project the zenith delays onto the slant paths and vice versa.

Current state of the art GNSS processing software use slightly more sophisticated tropospheric models where, in addition to the ZHD and ZWD parameters, horizontal gradients are included to describe a first order azimuthal variation in the tropospheric refractivity field above the antenna (Davis et al., 1993):

$$ \Delta L = \Delta L_{h}^{z} \times m_{h}(\epsilon) + \Delta L_{w}^{z} \times m_{w}(\epsilon) + (\vec{G} \cdot \vec{e}) \times m_{G}(\epsilon) \qquad (2.5) $$

where $\Delta L_h^z$ and $\Delta L_w^z$ are the zenith hydrostatic and wet delays, respectively, and $\vec{G}$ is the horizontal gradient vector, $\vec{e}$ signifies the unit vector pointing in the direction of the observed satellite, and $m_h$, $m_w$, and $m_G$ are the mapping functions for the ZHD, ZWD, horizontal gradient, respectively. Here $\epsilon$ represents the elevation angle.

In summary, the handling of tropospheric delay in GNSS data processing consists in mainly specifying: i) the mapping functions and the source of surface pressure or ZHD data (e.g. from an atmospheric model), ii) the rate at which ZWD and gradient parameters should be estimated (e.g. once per hour) to account for atmospheric variability, iii) the constraints on the temporal variation of the estimated parameters.
### 2.1.1.3 Conversion of tropospheric delay to IWV

The GNSS IWV product is derived in a post-processing step which consists in extracting the ZWD part from the Zenith Tropospheric Delay (ZTD) that is the standard output of the processing software. For this purpose, we don't use the ZHD that served for a priori correction in Equation (2.1) but a more accurate estimate if one is available. The reason is that with the current state of the art GNSS processing standards, the prior ZHD is taken from a coarse resolution atmospheric model (e.g. 2° × 2.5°, see next sub-section), which is also used to compute the mapping function models. A more accurate ZHD estimate can be computed from a high-resolution version of the atmospheric model. In this work we use ERA-Interim or ERA5, which have a grid resolution of 0.75° and 0.25°, respectively, in both dimensions. Once an accurate estimate of ZHD is subtracted from the ZTD, the ZWD component needs to be converted into IWV by application a "delay to mass" conversion factor, denoted $κ$. The definition of $κ$ stems from the ratio of IWV defined by Equation (1.3) and ZWD defined by:

$$ \Delta L_w^z = 10^{-6} R_v \int_0^\infty \rho_v(z) \times (k'_2 + \frac{k_3}{T_m}) dz \quad (2.6) $$

where $k'_2$ and $k_3$ are refractivity coefficients for the water molecule, and $R_v$ is the specific gas constant for water vapor. Various empirical values for the refractivity coefficients are available (Bevis et al., 1994). In this work we used the updated version proposed by Bock (2020a). Forming the ratio of Equations (1.3) and (2.6), we get:

$$ \kappa(T_m) = \frac{10^6}{R_v(k'_2 + \frac{k_3}{T_m})} $$

where $T_m$, referred to as the "weighted mean temperature", is given by:

$$ T_m = \frac{\int \rho_v(z) dz}{\int \frac{\rho_v(z)}{T(z)} dz} \quad (2.7) $$

The conversion of ZWD into IWV finally writes:

$$ IWV = \kappa(T_m) \times ZWD. \quad (2.8) $$

A rule of thumb for Equation (2.8) is $IWV = 155 \times ZWD$ when ZWD is expressed in m, or $IWV = ZWD/6.5$ when ZWD is expressed in mm, which means that 6.5 mm of wet delay at zenith converts to 1 kg m⁻² of integrated water vapour (Bock, 2020a). This rule is useful to quickly convert delay errors into IWV errors.

For a more accurate conversion of estimated ZWD data, knowledge of the $T_m$ parameter is necessary. Examination of (2.7) reveals that, although both the vertical profile of water vapour density, $ρ_v(z)$, and temperature, $T(z)$, appear in the definition of $T_m$, the former is both in the numerator and denominator of the fraction, and therefore the computed $T_m$ is only sensitive to second order to errors in $ρ_v(z)$. This suggests that $T_m$ could be parameterized as a function of temperature. Bevis et al. (1992) showed that $T_m$ can be predicted to a good accuracy from surface air temperature, and proposed an empirical model obtained by linear regression. Such a model is, however, highly site dependent, and needs to be recomputed for stations in different climatic regions. Moreover, the prediction based on surface temperature tend to map spurious diurnal variations into $T_m$. Because of these limitations, it is preferable and technically more convenient to use $T_m$ estimates based on a high resolution atmospheric model (Bock, 2020a).
In this work, the ZHD and $T_m$ conversion parameters were both computed from ERA-Interim or ERA5 vertical profiles at their nominal horizontal grid resolution. The uncertainty in these parameters is assumed to be below 2.5 mm in ZHD (equivalent to an error in surface pressure of 1 hPa) and below 3 K in $T_m$. The induced IWV equivalent uncertainty is 0.4 kg m⁻² and 0.3 kg m⁻², respectively, assuming a typical mid-latitude $T_m$ of 275 K and IWV of 25 kg m⁻².

### 2.1.1.4 Important data processing aspects

GNSS processing can be carried out using various software packages and processing options. In this subsection, we explore the significant factors influencing the estimated ZTD accuracy, such as processing mode, tropospheric modeling aspects, session duration, elevation cutoff angle, elevation-dependent observation weighting, and correction models for antenna phase center offsets and variations.

Scientific GNSS software packages are generally based on one of two processing modes: double difference (DD) positioning, which uses double-differenced phase observations from a network of stations, or precise point positioning, so-called PPP, which use the undifferenced phase observations (Teunissen and Montenbruck, 2017a, p. 13). DD processing is independent of external satellite clock products, while PPP relies on the precise satellite orbits and clocks. The increased availability of precise orbit and clock solutions in the late 1990s, thanks to the efforts of the IGS, has been instrumental. Moreover, PPP demands precise models to correct for systematic effects causing centimeter-level variations (e.g. the Higher-Order Ionospheric Delay Corrections, Site Displacement Effects, etc). An overview of the various model components and corrections in PPP applications is provided in Teunissen and Montenbruck (2017a, p. 727).

Several aspects of tropospheric modelling have been implemented differently depending on the software. This concerns especially the modelling of the temporal variations of the estimated ZWD and gradient parameters. For example, in the Bernese software, a piece-wise linear model is fitted using the least squares method (Dach et al., 2015). Conversely, the GIPSY software utilizes a random walk model (Zumberge et al., 1997), and the GAMIT software relies on the Gauss-Markov model (Herring et al., 2015). While GAMIT is also based on a least-squares estimation, the GIPSY software implements the Kalman Filtering technique. The latter approach makes it possible to estimate tropospheric parameters at each epoch of observations (e.g. every 5 min), but this requires to set adequate constraints in the random walk model.

Mapping functions are crucial elements of the tropospheric model as they represent the partial derivatives of the estimated parameters in Equation (2.5). Currently, the two widely used mapping functions are the first Vienna mapping function (VMF1) (Boehm et al., 2006b) and the Global Mapping Function (GMF) (Boehm et al., 2006a). They are both built on the concept of the Niell Mapping Function (NMF) (Niell, 1996) and the Isobaric mapping Function (IMF) (Niell, 2000) but use numerical weather model data which provide a more comprehensive description of the refractive index in space and time. These mapping functions are parameterized with a number of coefficients that need be extracted at the specific GNSS station locations by the processing software. In the case of the VMF1, the coefficients are provided with a 6-hourly time step from the operational ECMWF analyses on a 2° latitude by 2.5° longitude horizontal grid. This mapping function, along with its
more recent version VMF3, are the most accurate. Currently, VMF1 is recommended by the IERS Service.

Another important feature is the correction of the antenna phase-center offsets (PCOs) and variations (PCVs), both on the satellite and the receiving side (Teunissen and Montenbruck, 2017a, Chap. 17). The magnitude of the PCVs depends on the elevation angle of the satellite as observed from the ground, and it can vary between different types of antennas. These offsets and variations have a significant impact on the estimated ZTD delay. Use of a proper correction models is essential to address this issue. Earlier models (before 2006) were based on a relative calibration technique and introduced a systematic bias in the absolute ZTD estimates. The more recent, absolute calibrations, have resolved this problem and have been adopted since 2006 in the IGS. However, for a small number of old antennas, used in the IGS network and elsewhere, only relative calibration models are available. This issue is discussed in more detail in Section 2 of paper 1 for the subset of IGS stations used in our work.

ZTD and several other parameters estimated during the GNSS processing are sensitive to elevation dependent errors. Such errors may especially arise from errors in the mapping functions which represent a smooth relationship between zenith delay and slant delay, whereas the true atmosphere may be more complex in many situations (e.g. a passing meteorological frontal system). They may also be due to the application of a wrong or an imperfect antenna PCO/PCV model, as well as to multipath and signal interference (Teunissen and Montenbruck, 2017a, Chap. 15 and 16). Mapping function errors and multipath errors are typically larger at lower elevations. One commonly used strategy to minimize their impact is to increase the elevation cutoff angle or to down-weight low-elevation observations. However, low-elevation observations improve geometry and reduce formal ZTD errors. Thus, an optimum combination of elevation cutoff angle and weighting strategy needs to be adjusted to mitigate these errors.

The overall uncertainty in the ZTD and ZWD estimates resulting from the typical current state-of-the art GNSS processing strategies is estimated to 2-5 mm (Teunissen and Montenbruck, 2017a, Chap. 38). Using the rule of thumb introduced in previous sub-section, the ZTD uncertainty adds a further 0.3 to 0.8 kg m⁻² uncertainty to the derived IWV estimates.

### 2.1.1.5 GNSS tropospheric delay products used in this work

In this thesis, we used three distinct GNSS datasets, each representing a different generation of (re-)processing products: IGS repro1 (first generation), CODE REPRO2015 (second generation), and NGL repro3 (third generation).

The IGS repro1 dataset was produced between 2010 and 2011 by JPL/NASA in the framework of the first reprocessing campaign organized by IGS (Byun and Bar-Sever, 2009). It marks the first comprehensive reprocessing effort with collaborative contributions from multiple Analysis Centers. This reanalysis encompassed the entire history of GPS data collected by the IGS global network from January 1, 1995, to December 31, 2007, employing the best available models and methodologies at the time. JPL completed this dataset consistently by mid-2011, providing tropospheric products (ZTD and gradients) for 460 IGS stations, for the period 1995-2010 (inclusive). A quality-checked and IWV-converted sub-set of 120 stations with more than 15 years of measurements was prepared by Bock (2016) and used by Parracho (2017) and Quarello (2020).
The second dataset, CODE REPRO2015, was processed by the Center for Orbit Determination in Europe (CODE) in 2015, using the framework settled for second IGS reprocessing campaign (Dach et al., 2015). It includes 434 IGS stations and covers the period from January 1, 1994, to December 31, 2014. Until now, this dataset has been extended year after year with the operational CODE product (Dach et al., 2018), and quality-checked and IWV-converted estimates have been released by Bock (2019).

Details of the processing methodologies for the IGS repro1 and CODE REPRO2015 datasets can be found in section 2 of paper No. 1 where various important features described in the previous sub-section are summarized and discussed. To mention just one, we note that IGS repro1 was produced with GIPSY OASIS II software in PPP mode and CODE REPRO2015 with the Bernese GNSS software in a DD positioning mode. For the purpose of paper No. 1, we selected 81 common stations, ensuring that the time series in both datasets covered a minimum period of 15 years. The map in Figure 1 in this paper displays the available stations from both datasets and highlights the selected stations.

<table><caption>Table 2.1 – Summary of GNSS data processing parameters of the NGL repro3 dataset.</caption><tr><td>Software</td><td>GipsyX Version 1.0</td></tr><tr><td>Strategy</td><td>Precise Point Positioning (PPP)</td></tr><tr><td>Orbits, Clocks, ERPs</td><td>Daily Repro3.0 (JPL)</td></tr><tr><td>Reference Frame</td><td>IGS14</td></tr><tr><td>Antenna Calibration</td><td>igs14_www.atx</td></tr><tr><td>Window Length</td><td>30 hours</td></tr><tr><td>Elevation Cutoff Angle</td><td>7°</td></tr><tr><td>Observations</td><td>GPS</td></tr><tr><td>Observation Sampling Interval</td><td>5 minutes</td></tr><tr><td>Observation Weighting</td><td>σ² = 1/sin(e) where e is the elevation</td></tr><tr><td>Tropospheric model</td><td>ZHD and ZWD a priori:<br/>6-hourly ECMWF analysis (provided by TUV)<br/>VMF1 mapping functions (hydrostatic and wet).<br/>Random Walk model for ZWD and<br/>gradient parameters with constraints:<br/>3 mm h⁻¹/² (ZWD) and 0.3 mm h⁻¹/² (gradients)<br/>ZWD and gradient sampling: 5 min</td></tr><tr><td>Tropo files</td><td>ZTD and gradient estimates provided in<br/>SINEX files (0000, 0005, ...2345 UTC)</td></tr><tr><td>Coordinate estimates</td><td>Estimated once per 24 h</td></tr><tr><td>Ambiguity resolution</td><td>Fixed</td></tr></table>

The third dataset, NGL repro3, is sourced from the Nevada Geodetic Laboratory (NGL) (Blewitt et al., 2018). It represents the latest generation of GNSS products. NGL routinely collects and processes geodetic-quality
GPS observations from more than 20,000 stations worldwide, encompassing various regional and commercial networks, in addition to the widely used IGS network. The dataset covers the period from 1994 to present. Details of the processing relevant to ZTD estimation are provided in Table 2.1. One can note that NGL uses the GipsyX version 1.0 software from JPL in PPP mode together with JPL's repro3 final GPS orbits and clocks. The IWV data from this dataset were derived by applying a quality-checking and conversion procedure consistent with the other two datasets. The daily IWV estimates were used for the attribution task in the second paper (Chapter 4). In this work, the above-mentioned 81 stations from the CODE REPRO2015 dataset were selected as the main stations together with 628 nearby stations from the NGL dataset.

The NGL repro3 dataset uses a more modern version of JPL's GNSS processing software, and thus shares many similar features with the IGS repro1 dataset, including PPP mode, ZWD and gradient stochastic model, and elevation cutoff angles, among others. Its mapping functions are more modern, and in line with CODE REPRO2015. In contrast to IGS repro1, NGL repro3 adopted a strategy that down-weights low-elevation observations, making it less sensitive to multipath errors. The fixed ambiguity resolution approach is also an efficient way to reduce ZTD biases. NGL repro3 adopted also the more recent reference frame, IGS14 and the associated igs14.atx. The latter has resulted in a significant increase, up to 90%, in the number of stations within the IGS network that have now absolute calibrations (Rebischung and Schmid, 2016).

### 2.1.2 Reanalysis data

Climate reanalysis combines modern weather forecasting models with historical observations through data assimilation to deliver a global picture of weather and climate of the past as close to the reality as possible. This representation is known for its global coverage and temporal consistency, often described as "maps without gaps", extending from the Earth's surface to the top of the atmosphere (Fact Sheet: Reanalysis 2023). They utilize realistic models to extrapolate information from locally observed parameters to nearby locations and forward in time. Unlike operational analyses produced by Numerical Weather Prediction systems, reanalyses are created using a consistent assimilation system and the same numerical model, making them unaffected by changes in methods, model physics, or dynamics.

In our research, we focused on two specific reanalysis datasets produced by the European Centre for Medium-Range Weather Forecasts (ECMWF) for climate monitoring: ERA-Interim (Dee et al., 2011) and ERA5 (Hersbach et al., 2020). ERA-Interim is ECMWF's previous atmospheric reanalysis, based on a 2006 version of the Integrated Forecasting System (IFS) physical model. Its data span from 1979 until August 31, 2019, and no further updates are available. ERA5, on the other hand, represents the latest (Fifth) generation of ECMWF global reanalysis, based on a 2016 version of the ECMWF IFS. ERA5 offers data starting from 1959 continuously expanded in near real-time. ERA5 offers higher spatial and vertical resolutions (0.25° and 137 levels) compared to ERA-Interim (0.75° and 60 levels). Additionally, ERA5 benefits from increased satellite observations and improvements in its operational processes. As an example, ERA5 assimilated several satellite channels sensitive to humidity through an all-sky assimilation approach, in contrast to ERA-Interim, which employs a clear-sky assimilation approach. Moreover, ERA5 uses a 10-member ensemble 4D-Var system, providing estimates of uncertainties in the data, a feature not present in ERA-Interim (Hersbach, 2019). Various enhancements are expected to yield more accurate IWV estimates, as well as improvements in surface pressure
and air temperature, which are used as auxiliary data for IWV conversion.

We anticipate that the representativeness difference (Bock and Parracho, 2019) between GNSS and ERA5 will be smaller compared to that between GNSS and ERA-Interim, primarily due to ERA5's improved spatial resolution. This difference stems from the fact that the reanalysis grid-point values are representative of aerial averages while the GNSS observations are comparatively equivalent to point observations. It is important to note that we include ERA-Interim primarily because this work builds on the previous work by Parracho (2017) and Quarello (2020), and we wanted to investigate the impact of using different reanalysis data in the homogenization process.

It is essential to emphasize that the homogeneity of reanalysis data is uncertain. As mentioned earlier, ERA-Interim showed inhomogeneities when changes occurred in the observing system (Schroeder et al., 2016), and perhaps ERA5 as well (Allan et al., 2022). For these reason, we also consider the possibility of having change-points in the reanalysis used as a reference in forming the IWV difference series in the homogenization process.

## 2.2 Mathematical tools

### 2.2.1 GNSSseg segmentation method

In this section we recall briefly the segmentation method developed by Quarello et al. (2022) for the change-point detection problem in GNSS series (being the basis of the work of this thesis).

#### 2.2.1.1 Model

Quarello et al. (2022) modeled the series of IWV difference by a Gaussian independent random process $z = \{z_t\}_{t=1,...,n}$, such that

$$z_t \sim N(\mu_k + f_t, \sigma_{month}^2) \text{ if } t \in I_k^{mean} \cap I_{month}^{var} \text{ for } k = 1, ..., K,$$

where

* $z_t$ is the data at time $t$,

* K denotes the number of segments,

* $\mu_k$ represents the constant mean of the $k$th segment $I_k^{mean} = [[t_{k-1} + 1, t_k]]$, with $t_k$ indicating the position of the $k$th change-point with the convention $t_0 = 0$ and $t_K = n$,

* $f_t$ a function of time modeled by a Fourier series of order 4 accounting for the periodic bias:

$$f_t = \sum_{i=1}^{4} a_i \cos(2\pi t/T) + b_i \sin(2\pi t/T),$$

with $T = 365.25$ (days) and the coefficients $a_i$ and $b_i$ are weights for the series representing yearly, semiannual, ter-annual, and quarterly cycles, corresponding to $i$ ranging from 1 to 4 respectively,

* $\sigma_{month}^2$ stands for the monthly variance, consistent within the monthly interval $I_{month}^{var} = \{t, date(t) \in month\}$.
### 2.2.1.2 Estimation procedure

The proposed inference procedure is illustrated in Figure 2.2. This procedure consists of three main steps:

1. Estimation of the monthly variance $\sigma_{month}^2$ as proposed by Bock et al. (2019).

2. For a fixed number of segments $K$, the periodic bias function $f_t$ and both the change-points $t_k$ and the segment means $\mu_k$ for $k = 1, ..K$, are iteratively estimated by minimizing the residual sum of squares, denoted by $SSR_K$:

$$SSR_K = \sum_{k=1}^{K} \sum_{month} \sum_{t \in I_k^{mean} \cap I_{month}^{var}} \frac{(z_t - f_t - \mu_k)^2}{\hat{\sigma}_{month}^2}.$$

This iterative procedure stops when the relative change in the estimated values is smaller than a threshold.

3. This estimation procedure is performed for different values of $K$ (from 1 to $K_{max}$ that is the maximal considered number of segments) and thus returns $K_{max}$ different segmentations. The "best" one, i.e. the "best" number of segments is determined as follows:

$$\hat{K} = \operatorname{argmin}_K SSR_K + pen(K).$$

Different penalty functions (pen) have been considered and tested in Quarello et al. (2022): mBIC (Zhang and Siegmund, 2007), Lav (Lavielle, 2005), BM1 (Birgé and Massart, 2001) and BM2 (Lebarbier, 2005). The idea of adding a penalty to $SSR_K$ is to prevent overfitting, where the model attempts to have as many change-points as possible.
<figure><img src="image_8.png" /><figcaption>Figure 2.2 – Segmentation inference procedure (Quarello, 2020).</figcaption></figure>

### 2.2.1.3 Important features

Certain factors need to be considered when employing this method.

The foremost factor concerns the selection of the number of change-points, *K*. Specifically, the challenge lies in selecting the appropriate criterion. Quarello et al. (2022) tested the four different penalized criteria through both simplified simulations and real data analysis (IGS solution). The simulation findings revealed that Lav exhibited instability, showing a large variance in detection numbers, whereas other criteria were more conservative. When applied to the real data, mBIC tended to over-segment compared to the other criteria. This tendency was consistent when we applied this method to other solutions (CODE). To select the final change-points, the authors proposed a semi-automatic approach which combines a visual inspection and the results from criteria. This approach involved inspecting the monthly time series for each station, prioritizing
based on known changes in GNSS metadata and information from the TEQC software. Through this validation process, Quarello et al. (2022) concluded that BM1 had the highest acceptance rate at 57%, followed by Lav at 51% and BM2 at 45%. Consequently, in our study, we focus solely on the results derived from the BM1 criterion.

Second, an important consideration is that no dependence between the data is assumed in this segmentation model. However, we have observed a time-dependent aspect within our data, a phenomenon commonly encountered in climate data as well (see the subsection 4.1.3.3). Neglecting this feature could lead to an over-segmentation of the series, as discussed in Chakar et al. (2017).

We acknowledge the potential issue of confusion between the estimates of $f_t$ and $\mu_k$ particularly when dealing with seasonal variations. These variations can occasionally be mistakenly attributed to the mean rather than the periodic function, a situation that can arise with oversegmentation. Additionally, the absence of interannual variation within the model may lead to a false change in the mean to take it into account.

Finally, the method has been implemented in a R package named GNSSseg available on the CRAN, followed by a faster version called GNSSfast and available on github. In our study, segmentation results were obtained using the GNSSfast version. We adhered to all default settings, which include a maximum segment number of $K_{max} = 30$ and an iterative estimation threshold set at $10^{-1}$ (the stopping rule of the iterative estimation procedure), etc...

### 2.2.2 Linear Regression

In this thesis, we considered the linear regression model for various tasks, including trend estimation and testing the significance of the offset at a given change-point. In this section, we recall the two commonly used inference methods, namely, Ordinary Least Squares (OLS) and Generalized Least Squares (GLS) for this model.

The multiple linear model with $k$ explanatory numerical variables is given, for $n$ observations, by:

$$z_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \dots + \beta_k x_{ik} + \epsilon_i \quad \text{for } i = 1, \dots, n.$$

This model can be expressed in matrix form as:

$$z = X\beta + e,$$

where $X$ is the matrix of explanatory variables with size $[n \times (k + 1)]$, $y$ is the $[n \times 1]$ vector of responses, $e$ is the $[n \times 1]$ vector of errors and $\beta$ the $[(k + 1) \times 1]$ vector of coefficient regression, i.e. the unknown parameters:

$$X = \begin{bmatrix} 1 & x_{11} & x_{12} & \dots & x_{1k} \\ 1 & x_{21} & x_{22} & \dots & x_{2k} \\ \vdots & \vdots & \vdots & \dots & \vdots \\ 1 & x_{n1} & x_{n2} & \dots & x_{nk} \end{bmatrix}, y = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}, e = \begin{bmatrix} \epsilon_1 \\ \vdots \\ \epsilon_{n} \end{bmatrix}, \beta = [\beta_0 \ \beta_1 \ \beta_2 \ \dots \ \beta_k] \text{ respectively.}$$

The least-squares method is the most popular inference method for linear regression. It consists in finding the
value of $\beta$ that minimizes the sum of squared errors given by

$$S(\beta) = (z - X\beta)^T (z - X\beta).$$

The solution, denoted by $\hat{\beta}$, is obtained by setting the (matricial) derivative of $S(\beta)$ with respect to $\beta$ equal to 0. This leads to the equation:

$$X^T X \hat{\beta} = X^T z.$$

When $X^T X$ is invertible, we obtain the following least-square (LS) estimator

$$\hat{\beta} = (X^T X)^{-1} X^T z. \qquad (2.9)$$

To study the properties of this estimator, let us substitute $z$ into the previous equation. Hence:

$$\begin{aligned} \hat{\beta} &= (X^T X)^{-1} X^T (X\beta + e) \\ &= (X^T X)^{-1} X^T X\beta + (X^T X)^{-1} X^T e \\ &= \beta + (X^T X)^{-1} X^T e. \end{aligned}$$

We can deduce some properties of the LS estimator. First, the LS estimator is unbiased as long as the noise is centered. Second, the uncertainty of the LS estimator depends on the variance of the noise (or equivalently of the variance of the response):

$$var(\hat{\beta}) = (X^T X)^{-1} (X^T var(z) X) (X^T X)^{-1}. \qquad (2.10)$$

Assuming that $e \sim N(0, \sigma_0^2 I_n)$, where $I_n$ is an identity matrix with size $n$, or at least that the errors are uncorrelated and equally variable, leads to the LS estimator with a reduced variance

$$var(\hat{\beta}) = \sigma_0^2 (X^T X)^{-1}.$$

This estimator is known as the Best Linear Unbiased Estimator (BLUE) of $\beta$ (Gauss-Markov theorem). BLUE implies that the variance of the estimator is the smallest compared to other unbiased linear estimators. This solution is the familiar Ordinary Least Squares (OLS), often denoted $\hat{\beta}^{\text{OLS}}$.

However, the OLS is no longer the BLUE when one of the two assumptions (independence and homoscedasticity) is violated, i.e. $var(e) = var(z) = \Sigma_0 \neq \sigma_0^2 I_n$ where $\Sigma_0$ is a positive semidefinite matrix. Indeed, in this case, we get

$$E(\hat{\beta}^{\text{OLS}}) = \beta,$$

the OLS is still unbiased, but according to equation (2.10), we get

$$var(\hat{\beta}^{\text{OLS}}) = (X^T X)^{-1} (X^T \Sigma_0 X) (X^T X)^{-1} \neq \sigma_0^2 (X^T X)^{-1}. \qquad (2.11)$$

A serious consequence is that statistical inferences (as tests) based on the standard OLS estimation results
become invalid.

The adapted framework in this case is the General Least Squares (GLS) method which aims to minimize the following sum of squares:

$$ (z - X\beta)^T \Sigma_0^{-1} (z - X\beta). $$

The solution of this minimization is called the GLS estimator given by

$$ \hat{\beta}^{\text{GLS}} = (X^T \Sigma_0^{-1} X)^{-1} X^T \Sigma_0^{-1} z, \qquad (2.12) $$

with variance

$$ \mathrm{var}(\hat{\beta}^{\text{GLS}}) = (X^T \Sigma_0^{-1} X)^{-1}. $$

Note that the GLS method is equivalent to the OLS method applied on the transformed model:

$$ G z = GX\beta + Ge, $$

with $G = \Sigma_0^{-1/2}$. The GLS estimator is the BLUE of $\beta$.

In practice, $\Sigma_0$ is unknown. Depending of the study objective, various method can be used. Two objectives can be considered: (i) when we are interesting in the estimator itself, the popular method is the feasible generalized least-squares (FGLS) estimator which consists in substituting a good estimator of $\Sigma_0$ in the (2.12) and (ii) when we are rather interesting in statistical inferences (as tests), the OLS estimator can be considered with a correction of its variance. More precisely, we can estimate the term $(X^T \Sigma_0 X)$ in the variance (see equation (2.11)) using a robust covariance (HAC) estimator.

### 2.2.3 Linear stochastic models

#### 2.2.3.1 Definition

In modeling time series, specific stationary stochastic processes that are valuable include the autoregressive (AR), moving average (MA), and mixed autoregressive-moving average (ARMA) processes (Box et al., 2016; Shumway and Stoffer, 2017). A generalization of the latter to autoregressive integrated moving-average (ARIMA) models provides also a useful class of nonstationary processes. While nonstationary processes are frequent in many application fields, in this work we are mainly concerned with stationary processes.

The autoregressive model of order $p$, or more succinctly, the AR($p$) process, is defined by the following equation Box et al. (2016, p. 52):

$$ z_t = \phi_1 z_{t-1} + \dots + \phi_p z_{t-p} + a_t, \qquad (2.13) $$

where $z_t$, assumed to be stationary and of zero mean, denotes the dependent variable at time $t$, $z_{t-1}$ to $z_{t-p}$ represent its previous values up to lag $p$, and $\phi_1$ to $\phi_p$ are the set of model parameter. The term $a_t$ represents an independent random sequence, or white noise, also known as the innovation. The latter follows a normal distribution with mean $\mu_a = 0$ and variance $\sigma_a^2$, i.e. $a_t \sim N(0, \sigma_a^2)$. To account for a non-zero mean, $\mu$, in $z_t$, a constant can be added to the right-hand side of equation (2.13), or $\mu$ can be subtracted from $z_t$ and all its lagged
terms in the equation.

The moving average process of order $q$, or MA($q$), is defined by the following equation Shumway and Stoffer (2011, p. 90):

$$z_t = a_t + \theta_1 a_{t-1} + \dots + \theta_q a_{t-q}, \quad (2.14)$$

where $a_{t-1}$ to $a_{t-q}$ represent the lagged noise terms up to lag $q$, and $\theta_1$ to $\theta_q$ correspond to the model coefficients associated with these lagged error terms. Note that in equation (2.14), we adopted a positive sign in the sum of the right-hand side terms, following the notation of Shumway and Stoffer (2011) and the convention in R, which is opposite to the one used in Box et al. (2016).

Finally, the definition of the ARMA($p, q$) process, as formulated by Shumway and Stoffer (2011, p. 92)), is:

$$z_t = \phi_1 z_{t-1} + \dots + \phi_p z_{t-p} + a_t + \theta_1 a_{t-1} + \theta_q a_{t-q}. \quad (2.15)$$

This expression represents a combination of the autoregressive model of order $p$ and the moving average model of order $q$.

The AR($p$) model implies that observations at time $t$ are influenced by $p$ previous observations with weights $\phi_1, \phi_2, \dots, \phi_p$, in addition to an independent white noise term at time $t$. On the other hand, the MA($q$) model suggests that observations at time $t$ are affected by $q$ previous noise terms with weights $\theta_1, \theta_2, \dots, \theta_q$.

To express the time series lags in a concise mathematical form, the backward shift operator $B$ is utilized, which is defined as follows:

$$Bz_t = z_{t-1}.$$

The AR($p$) model can thus be reformulated as:

$$(1 - \phi_1 B - \dots - \phi_p B^p) z_t = a_t.$$

Introducing the autoregressive polynomial $\phi(B) = (1 - \phi_1 B - \dots - \phi_p B^p)$, we have:

$$\phi(B)z_t = a_t.$$

Similarly, the MA($q$) model can be rewritten as:

$$z_t = \theta(B)a_t.$$

where the moving average polynomial is denoted as $\theta(B) = (1+\theta_1 B+\dots+\theta_q B^q)$. Finally, the mixed ARMA($p, q$) model is represented as:

$$\phi(B)z_t = \theta(B)a_t.$$
### 2.2.3.2 Important features

A central feature in the development of time series models is the assumption of stationarity. Stationarity implies that the mean, $E[z_t]$, variance, $var[z_t] = E[(z_t - \mu)^2]$, and autocovariance between time $t$ and $t-h$, $cov[z_t, z_{t-h}] = E[(z_t - \mu)(z_{t-h} - \mu)]$, are the same for all times $t$ Shumway and Stoffer (2011, p. 23). For a stationary series, we have thus:

$$E[z_t] = \mu, \quad var[z_t] = \sigma_z^2, \quad \text{and} \quad cov[z_t, z_{t-h}] = \gamma(h).$$

Note that the autocovariance depends only on the time lag $h$.

While the MA($q$) model is always stationary, the AR($p$) model does not always exhibit stationarity. It is shown in Box et al. (2016, pp. 54–55) that the AR($p$) model achieves stationarity only when the roots of the autoregressive polynomial lie outside the unit circle.

Another crucial property for a linear process is invertibility. It relates to the capability of representing the error term $a_t$ as a function of $z_t$, making it algebraically equivalent to a converging infinite order AR model. The AR($p$) model is generally invertible, while the MA($q$) model is not always so. The invertibility of the MA($q$) model is achieved on the condition that $\theta(B)^{-1} \neq 0$, which implies that the roots of the moving average polynomial lie outside the unit circle (Box et al., 2016, pp. 68–69). For an ARMA($p$, $q$) process to be stationary and invertible, both conditions must be satisfied simultaneously.

To illustrate these concepts, the following sub-section will describe the properties of three fundamental models, the AR(1), MA(1), and ARMA(1,1), that will be extensively used in this work. For simplicity, we will henceforth denote $\phi_1$ and $\theta_1$ as $\phi$ and $\theta$, respectively, since the models are of first order.

### 2.2.3.3 Fundamental models

#### AR(1) model

The zero-mean first-order autoregressive model is expressed as:

$$z_t = \phi z_{t-1} + a_t,$$

where the stationarity condition restricts $|\phi| < 1$ and where $\{a_t\}_t$ i.i.d. $\sim N(0, \sigma_a^2)$. When $|\phi| = 1$, it transforms into the (non-stationary) random walk model.

We can easily verify that the mean is zero by writing:

$$E[z_t] = E[\phi z_{t-1}] + E[a_t] = \phi E[z_{t-1}] = 0,$$

given that $E[a_t] = 0$ and that, due to stationarity, the mean of the process must be constant, hence $E[z_t] = E[z_{t-1}]$.
The variance can be expressed as:

$$var[z_t] = E[z_t^2] = \frac{\sigma_a^2}{1 - \phi^2},$$

given that due to stationarity, the variance must be constant, i.e. $var[z_t] = \phi^2var[z_{t-1}] + \sigma_a^2 = var[z_{t-1}]$.

The autocovariance function (ACVF) at lag $h$ writes (Shumway and Stoffer, 2011, p. 86):

$$\gamma(h) = \operatorname{cov}[z_t, z_{t-h}] = \frac{\sigma_a^2 \phi^h}{1 - \phi^2}, h \ge 0,$$

with $\gamma(h) = \gamma(-h)$, and the Autocorrelation function (ACF) at lag $h$ is:

$$\rho(h) = \operatorname{corr}[z_t, z_{t-h}] = \frac{\gamma(h)}{\gamma(0)} = \phi^h, h \ge 0,$$

with the property $\rho(h) = \rho(-h)$.

The ACF tails off exponentially and rapidly approaches 0 due to the condition $|\phi| < 1$. When $\phi = 0$, the process becomes white noise, whose ACF is $\rho(h) = 0$ for all lags except for lag $h = 0$ where its value is $\rho(0) = 1$.

An important property of linear processes is causality, the property by which $z_t$ depends only on past values, $z_{t-1}$, $z_{t-2}$, $a_t$, $a_{t-1}$... Time series issued from the physical measurements we are concerned with in this study are always causal.

In the case of the AR(1) process, causality is equivalent to, or follows from, stationarity. This can be seen by arranging the recursive equation (2.13) as a function of $a_t$, $a_{t-1}$..., or, equivalently, by inverting the polynomial equation, $(1 - \phi B)z_t = a_t$, as $z_t = (1 - \phi B)^{-1}a_t = \sum_{j=0}^{\infty} \phi^j a_{t-j}$. This alternative representation of the AR(1) process is known as its "infinite MA representation", which is also basically the form of any causal linear process. Stationarity of this process requires that the infinite sum is converging, which is obtained when $|\phi| < 1$.

### MA(1) model

The zero-mean first-order moving average model is represented as:

$$z_t = \theta a_{t-1} + a_t,$$

where $\{a_t\}_t$ i.i.d. ~ $N(0, \sigma_a^2)$. Again, we can check that the mean of this process is zero:

$$E[z_t] = \theta E[a_{t-1}] + E[a_t] = 0, \qquad (2.16)$$

and the variance is

$$var[z_t] = \theta^2 var[a_{t-1}] + var[a_t] = (1 + \theta^2)\sigma_a^2. \qquad (2.17)$$

by virtue of the independence and stationarity of the white noise.
The ACVF at lag $h$ writes (Shumway and Stoffer, 2011, p. 90):

$$ \gamma(h) = \begin{cases} (1 + \theta^2)\sigma_a^2, & h = 0, \\ \theta^2\sigma_a^2, & h = 1, \\ 0, & h > 1. \end{cases} \qquad (2.18) $$

and the ACF is:

$$ \rho(h) = \begin{cases} 1, & h = 0, \\ \frac{\theta}{1+\theta^2}, & h = 1, \\ 0, & h > 1. \end{cases} $$

The ACVF and ACF of the MA(1) process are equal to zero for all lags higher than 1. For the general MA($q$), this holds for all lags larger than $q$. The MA($q$) process is thus always stationary.

One caveat of the MA($q$) model is the non-uniqueness of its ACVF and ACF. For the MA(1) model, for example, $\rho(h)$ is the same for $\theta$ and $1/\theta$, and $\gamma(h)$ is the same for the pair $(\theta, \sigma_a^2 = 1)$, and $(1/\theta, \sigma_a^2 = 1/\theta^2)$. In practice, one chooses the model which satisfies the invertibility condition, $|\theta| < 1$ (Shumway and Stoffer, 2011, p. 92).

In the backward operator notation, the MA(1) model can be expressed as $z_t = (1 + \theta B)a_t$. When $|\theta| < 1$, the inverse of the operator $(1 + \theta B)$ is given by:

$$ (1 + \theta B)^{-1} = \sum_{j=0}^{\infty} (-\theta)^j B^j. $$

This infinite sum converges for all $|B| \le 1$. Consequently, the MA(1) model can be equivalently written as:

$$ (\sum_{j=0}^{\infty} (-\theta)^j B^j) z_t = a_t. $$

This is the infinite AR representation of the MA(1) model, and it only exists when $|\theta| < 1$, i.e. when the model is invertible.

### ARMA(1,1) model

The zero-mean first-order autoregressive moving average model is given by:

$$ z_t = \phi z_{t-1} + \theta a_{t-1} + a_t, $$

where $\{a_t\}_t$, i.i.d. ~ N(0, $\sigma_a^2$). To be causal and unique, it necessitates both the stationarity condition of AR(1) and the invertibility condition of MA(1). Under these conditions, the ACVF at lag $h$ writes (Shumway and Stoffer, 2017, p. 92):

$$ \gamma(h) = \begin{cases} \frac{\sigma_a^2 (1 + 2\theta\phi + \theta^2)}{1 - \phi^2}, & h = 0, \\ \frac{\sigma_a^2 (1 + \theta\phi)(\theta + \phi)\phi^{h-1}}{1 - \phi^2}, & h \ge 1. \end{cases} $$
The ACF is:

$$ \rho(h) = \begin{cases} 1, & h = 0, \\ \frac{(1+\theta\phi)(\theta+\phi)}{1+2\theta\phi+\theta^2}\phi^{h-1}, & h \ge 1. \end{cases} $$

The general pattern of the ACF of the ARMA(1,1) model closely resembles that of the AR(1) model, both exhibiting exponential decay with h. Because of this similarity, relying solely on the ACF proves inadequate for distinguishing between the two models. This consideration will lead to the use of the Partial autocorrelation function (PACF).

Due to the properties of the AR(1) and MA(1), the ARMA(1,1) process can be represented either as an infinite MA process or as an infinite AR process, and these properties can actually be generalized to the ARMA(p, q) model as well.

A prevalent issue of the ARMA(p, q) models is the parameter redundancy, or over-parameterization. In the case of the ARMA(1,1) model, this can be evidenced by reformulating in the backward operator form:

$$ (1 - \phi B)z_t = (1 + \theta B)a_t. $$

It is obvious that in the scenario where $\phi = -\theta$, $z_t = a_t$ which is actually a white noise process. As this example points out, we might fit an ARMA(1, 1) model to white noise data and find that the parameter estimates are significant. If we were unaware of parameter redundancy, we might claim the data are correlated when in fact they are not (Shumway and Stoffer, 2017, p. 84). This situation highlights a challenge involving the potential confusion between the white noise model and the ARMA(1,1) model that may arise during the process of model identification. This issue can be solved, in the general ARMA(p, q) model, by reducing the common factors in the $\phi(B)$ and $\theta(B)$ polynomials (Shumway and Stoffer, 2017, p. 85).

### 2.2.3.4 Note on the PACF

The Partial Autocorrelation Function (PACF) is a tool that helps to identify the order of an AR(p) model and, in combination with the ACF, helps to distinguish between AR, MA, and ARMA models. The idea is to construct a function that cuts off after lag p for an AR(p) model, like ACF cuts off after lag q for a MA(q) model. The principle is the following (Shumway and Stoffer, 2017, p. 96): if X, Y, and Z are random variables, then the partial correlation between X and Y given Z is obtained by regressing X on Z to obtain $\hat{X}$, regressing Y on Z to obtain $\hat{Y}$, and then calculating $\rho_{XY|Z} = corr(X - \hat{X}, Y - \hat{Y})$. The PACF of a stationary process, $z_t$, denoted $\phi_{hh}$, for $h = 1, 2, \dots$, is (Shumway and Stoffer, 2017, p. 97):

$$ \begin{cases} \phi_{11} = corr(z_{t+1}, z_t) = \rho(1); \\ \phi_{hh} = corr(z_{t+h} - \hat{z}_{t+h}, z_t - \hat{z}_t), & h \ge 2, \end{cases} $$

where $\hat{z}_{t+h}$ is the regression of $z_{t+h}$ on $z_{t+h-1}, z_{t+h-2}, \dots, z_{t+1}$.
In the case of a stationary AR(1) process of parameter $\phi$, we have:

$$ 
\begin{cases}
\phi_{11} = \rho(1) = \phi, \\
\phi_{22} = \operatorname{corr}(z_{t+2} - \hat{z}_{t+2}, z_t - \hat{z}_t) = \operatorname{corr}(a_{t+2}, z_t - \hat{z}_t) = 0, \\
\phi_{hh} = 0, \quad h > 2/
\end{cases}
$$

This indicates that for an AR(1) model, the PACF cuts off after the first lag. The property is the same for the AR($p$) models with the PACF cutting off after lag $p$.

For an invertible MA(1) model of parameter $\theta$, the PACF writes (Shumway and Stoffer, 2011, p. 107):

$$ 
\begin{cases}
\phi_{11} = \frac{\theta}{1+\theta^2}, \\
\phi_{22} = -\frac{\theta^2}{1+\theta^2+\theta^4}, \\
\phi_{hh} = -\frac{(-\theta)^h(1-\theta^2)}{1-\theta^2(h+1)}, \quad h > 2.
\end{cases}
$$

From this equation, it is clear that the PACF for the MA(1) model will not cut off, contrary to its ACF, but instead it will tail off like the ACF of an AR(1) model. This behavior can be understood by recognizing that an invertible MA process can be represented as an infinite AR process, preventing the PACF from cutting off at any specific lag. The PACF of a causal and invertible ARMA(1,1) model exhibits the tailing off like the MA(1) model.

<figure><img src="image_9.png" /><figcaption>Figure 2.3 – Theoretical ACF and PACF of AR(1) model with $\phi = 0.3$, MA(1) model with $\theta = 0.3$, and ARMA(1,1) model with (a) $\phi = 0.6$ and $\theta = -0.3$, and (b) $\phi = 0.3$ and $\theta = 0.2$.</figcaption></figure>

Figure 2.3 presents the ACF and PACF for the three model types with typical $\phi$ and $\theta$ values encountered in the real data that are discussed in this work (see Chapter 4). The three model types can be clearly distinguished. For the AR(1) model, the ACF tails off, while the PACF cuts off after lag $p = 1$. For the MA(1) model, the ACF cuts off after lag $q = 1$, while the PACF tails off. And for the two ARMA(1,1) models, both the ACF and PACF
tail off. However, depending on the sign of $\phi + \theta$, the PACF oscillates around zero or not. All these general behaviours are also valid for higher order models, AR(p), MA(q), and ARMA(p, q).

### 2.2.3.5 Note on spectral properties of stationary models

Time series often exhibit periodic variations which suggest that they can be approximated by a limited-order Fourier series. The concept can be generalized to non-periodic time series and linear stochastic processes by introducing a continuous distribution of frequencies (instead of harmonics considered in the simple Fourier Series). A fundamental tool in spectral analysis of stationary processes is the spectral density, also called power spectral density (PSD), $f(\omega)$, which is defined as the Fourier Transform of the ACVF (Shumway and Stoffer, 2017, p. 173):

$$f(\omega) = \sum_{h=-\infty}^{\infty} \gamma(h)e^{-2\pi i \omega h}, -1/2 \le \omega \le 1/2, \quad (2.19)$$

where $\omega$ is the continuous frequency variable (in cycles per point). Because $z_t$ is stationary, the inverse transform of $f(\omega)$ exists, and yields the ACVF:

$$\gamma(h) = \int_{-1/2}^{1/2} f(\omega)e^{2\pi i \omega h} d\omega, h = 0, \pm 1, \pm 2, \dots$$

i.e. the PSD and the ACVF are Fourier Transforms pairs.

The fact that $\gamma(h)$ is non-negative ensures that $f(\omega) \ge 0$ for all $\omega$, and it follows also from equation (2.19) that $f(\omega) = f(-\omega)$, and that $\gamma(0) = \text{var}[z_t] = \int_{-1/2}^{1/2} f(\omega)d\omega$.

The PSD of some special cases of interest are given below (Box et al., 2016): for $-1/2 \le \omega \le 1/2$,

- White Noise (WN): $f(\omega) = \sigma_a^2$,

- AR(1): $f(\omega) = \frac{\sigma_a^2}{1+\phi^2-2\phi\cos(2\pi\omega)}$,

- MA(1): $f(\omega) = \sigma_a^2 \times [1 + \theta^2 + 2\theta\cos(2\pi\omega)]$,

- ARMA(1,1): $f(\omega) = \sigma_a^2 \times \frac{(1+\theta^2+2\theta\cos(2\pi\omega))}{1+\phi^2-2\phi\cos(2\pi\omega)}$.

Figure 2.4 displays the PSDs for the same four models previously depicted in Figure 2.3. While AR(1) and MA(1) share similarities such as higher power in the lower frequencies and finite power at the limits $\omega = 0$ and $\omega = 1/2$. Note that the negative frequency range is usually not plotted due to the symmetry of the PSD ($f(\omega) = f(-\omega)$). The main difference between the AR(1) and MA(1) plots is that the PSD of the MA(1) is flatter at low frequencies and the PSD of the AR(1) decreases faster in the medium frequency range. Comparatively, the PSD of the first ARMA(1,1) model starts higher than the AR(1) and MA(1) and decreases faster the former two and than the second ARMA model. These distinguishable features can help to identify the different models from a visual inspection, as an alternative or complementary way to the ACF and PACF discussed above.

In practice, the difficulty is that real data do not necessarily follow a simple stochastic model and the usual ACF, PACF, and PSD estimators are subject to errors, making the identification more difficult.
<figure><img src="image_10.png" /><figcaption>Figure 2.4 – Theoretical spectral density of AR(1) model with $\phi = 0.3$, MA(1) model with $\theta = 0.3$, and ARMA(1,1) model with (a) $\phi = 0.6$ and $\theta = -0.3$, and (b) $\phi = 0.3$ and $\theta = 0.2$. In all cases, the variance of the innovation is $\sigma_a^2 = 1$.</figcaption></figure>

### 2.2.3.6 Combination of models

In practice, we may often consider that an observed time series is the sum of two or more independent series, e.g. a physical signal following an AR(1) process summed with an additive white noise (WN). The properties of the sum of ARMA-type models was first studied by Granger and Morris (1976). They showed that, if $X_t \sim \text{ARMA}(p, m)$, $Y_t \sim \text{ARMA}(q, n)$, and $Z_t = X_t + Y_t$, then $Z_t \sim \text{ARMA}(x, y)$, where $x \le p + q$ and $y \le \max(p+n, q+m)$. The result can be generalized to the sum of N ARMA models:

$$ \sum_{j=1}^{N} \text{ARMA}(p_j, m_j) = \text{ARMA}(x, y), $$

where

$$ x \le \sum_{j=1}^{N} p_j, $$

and

$$ y \le \max(x - p_j + m_j, j = 1, \dots, N). $$

Some specific combinations are of particular interest, namely:

1. $\text{AR}(1) + \text{WN} = \text{ARMA}(1,1)$,

2. $\text{AR}(1) + \text{AR}(1) = \text{ARMA}(2,1)$,

3. $\text{MA}(1) + \text{WN} = \text{MA}(1)$,

4. $\text{MA}(1) + \text{MA}(1) = \text{MA}(1)$,

5. $\text{ARMA}(1,1) + \text{WN} = \text{ARMA}(1,1)$,

6. $\text{AR}(1) + \text{MA}(1) = \text{ARMA}(1,2)$.

Hence, it is quite likely that real data resulting from a combination of processes give rise to ARMA models.

Granger and Morris (1976) also address the question whether a given specific ARMA($p, q$) model could have arisen from simpler models. Of interest to us is, especially, the case whether an ARMA(1,1) results from the
combination of an AR(1) with WN. They show that the answer is positive when a realizability condition is met with the ARMA(1,1) model parameters, $\phi$ and $\theta$, namely, if:

$$ \frac{1}{1 - \phi^2} > \frac{\rho_1}{\phi} \geq 0, $$

where $\rho_1 = \frac{\theta}{1+\theta^2}$ is the ACF at lag $h = 1$ of the MA(1) part of the ARMA(1,1) process, then the ARMA(1,1) = AR(1) + WN.

### 2.2.3.7 Model identification

The model identification for general ARMA($p, q$) models aims at determining the orders $p$ and $q$. The widely used Box-Jenkins method relies extensively on a visual analysis of ACF and PACF plots (Box et al., 2016, pp. 180–185). PSD plots are sometimes used as well, although these are more commonly used for building transfer function and multivariate models (Box et al., 2016, Chap. 12).

In recent years, endeavors have been made to automate identification processes through estimation-based methods (Hyndman and Khandakar, 2008; Koreisha and Pukkila, 1995). The underlying concept involves fitting a range of potential ARMA models (different values of couple $(p, q)$) and then selecting the optimal one using a penalized version of the fitted criterion. The most common penalty criteria are the Akaike information criterion (AIC), bias corrected AIC (AICc), and Bayesian information criterion (BIC). The AIC is defined for a model $m$ as:

$$ AIC(m) = -2 \ln(\hat{L}(m)) + 2 D(m), $$

and the BIC as:

$$ BIC(m) = -2 \ln(\hat{L}(m)) + D(m) \ln(n), $$

where $\hat{L}(m)$ denotes the estimated likelihood of model $m$, $D(m)$ is its number of parameters to be estimated, and $n$ is the sample size. Since the BIC penalty is higher than the AIC one, BIC tends to select simplest or parsimonious models compared to AIC (i.e. small values of $p$ and $q$). Note that for ARMA framework, AIC is known to overestimate the orders (Shibata, 1976).

In our study, we consider the BIC criterion and use the "auto.arima" function with $d = 0$ of the R package "forecast" (Hyndman et al., 2018).

### 2.2.3.8 Parameter estimation

To estimate the coefficient parameters in ARMA models, the well known methods proposed in the literature are the maximum likelihood, the conditional or unconditional ones (or equivalently the conditional or unconditional Sum of Squares), and for the specific case of AR models, the method of moments using the Yule-Walker equations. The exact or unconditional likelihood leads to a complex optimization problem since it is non-linear according to the parameters and thus requires the use of iterative algorithms. The conditional likelihood (conditionally to the initial observations) method is a way to simplify the problem. However, even if for very large samples, the both likelihood estimators are equivalent with the same asymptotic distribution (Hamilton 1994, p. 126), the conditional likelihood can result in biased estimates for relatively short series (Box et al., 2016, p. 526). Because of this, the use of the unconditional likelihood function is typically recommended for models
with moving average terms. For a complete presentation and discussion of these inference methods, see Box et al. (2016, Chap. 7) or Brockwell et al. (1991, Chap. 8). In our study, we consider the unconditional inference likelihood approach and use the "arima" function in R.

The large sample variance for the maximum likelihood estimates of an AR(1), MA(1) and ARMA(1,1) models are (Box et al., 2016, pp. 233–238):

AR(1) model:

$$var[\hat{\phi}] \approx \frac{1 - \phi^2}{n}, \qquad (2.20)$$

MA(1) model:

$$var[\hat{\theta}] \approx \frac{1 - \theta^2}{n}, \qquad (2.21)$$

ARMA(1,1) model:

$$var[\hat{\phi}] \approx \frac{(1 + \phi\theta)^2(1 - \phi^2)}{n(\phi + \theta)^2}, \qquad (2.22)$$

$$var[\hat{\theta}] \approx \frac{(1 + \phi\theta)^2(1 - \theta^2)}{n(\phi + \theta)^2}. \qquad (2.23)$$

These equations underscore the relationship between the variance of coefficient estimates, the length of series, $n$, and the values of the coefficient, $\phi$ and $\theta$. As one could guess, the variances are inversely proportional to $n$ and are decreasing when the parameters $\phi$ or $\theta$ increase.

To quantify, let's consider an AR(1) process with a coefficient value of $\phi = 0.3$ and a series of length $n = 1000$. According to equation (2.20), $var[\hat{\phi}] = (0.03)^2$, signifying a relative uncertainty in $\hat{\phi}$ of 0.03/0.3=0.1, i.e. 10%. However, if we maintain the same series length while increasing the coefficient value to $\phi = 0.6$, the variance reduces to $(0.025)^2$, which corresponds to a relative uncertainty of 4%.

Now let's consider a scenario where the model is identified as an ARMA(1,1) but the true model is AR(1), with $\phi=0.3$. In this case, $var[\hat{\phi}]$, as given by equation (2.22), amounts to $var[\hat{\phi}] = (0.042)^2$, which is 1.96 times larger than the variance of the AR(1) model of $(0.03)^2$. This illustrates the over-fitting caveat and emphasizes the necessity of identifying the appropriate stochastic model rather than relying on a more general model.

### 2.2.4 Classification

Classification is a supervised machine learning method whose aim is to predict the membership of an individual to a given class of the population, summarized by the class number $y$ called the label, from the information available on this individual $x$. Solving a classification problem amounts to construct a classifier, i.e. a function $\psi$ such that $\psi(x)$ represents the prediction of $y$ given $x$. We naturally want to build a high-performance classifier, whose prediction error rate is as low as possible. Assume that $(x, y)$ is a realization of the random vector $(X, Y)$ of (unknown) distribution $P$. The classifier makes an error when $\psi(X) \neq Y$, and the quality of the classifier $\psi$ can be measured thus by the probability of misclassification:

$$\mathcal{R}_P(\psi) = E[\mathbb{1}_{\psi(X) \neq Y}] = P(\psi(X) \neq Y),$$
This probability is also known as the risk of $\psi$. The ideal or optimal classifier will be the one that minimizes the risk $\mathcal{R}_P(\psi)$ denoted $\psi^*$. This classifier is particular and is called the Bayes classifier ($\psi^* = \arg\max_k P(Y = k|X = x)$). However, $\psi^*$ depends on the distribution $P$ of $(X, Y)$ which is unknown in practice. The objective of the supervised classification is thus to construct a learning rule or classifier from a sample $D_n = \{(X_1, Y_1), \dots, (X_n, Y_n)\}$ with the same distribution as $(X, Y)$ and with performance close to $\psi^*$ in terms of risk. This classifier is denoted $\hat{\psi}(D_n)$.

There are four popular learning algorithms in the literature: the Classification and Regression Tree (CART), the Random Forest (RF), the k nearest neighbors (k-NN) and the Linear Discriminant Analysis (LDA). The first three algorithms are non-parametric, in contrast to LDA. We describe these algorithms in subsection 2.2.4.1. The performance of a learning algorithm is summarized by the risk of its resulting classifier, $\mathcal{R}_P(\hat{\psi}(D_n))$. However, again the distribution $P$ being unknown, this risk can not be calculated. We can estimate it by its empirical version on the sample $D_n$:

$$ \hat{\mathcal{R}}_n(\psi; D_n) = \frac{1}{n} \sum_{i=1}^{n} \mathbb{1}_{\{\psi(X_i) \neq Y_i\}}, \qquad (2.24) $$

This is referred to as the empirical risk, representing the number of errors the classifier $\psi$ makes on the sample $D_n$. It serves as an unbiased estimator of the risk associated with $\psi$, i.e. $E[\hat{\mathcal{R}}_n(\psi; D_n)] = \mathcal{R}_P(\psi)$.

The quality of the classifier $\hat{\psi}(D_n)$ is thus given by $\hat{\mathcal{R}}_n(\hat{\psi}(D_n); D_n)$. However, this measure of quality is biased because it is associated to an optimistic estimation of the risk $\mathcal{R}_P(\hat{\psi}(D_n))$. As clearly indicated by these notations, the sample $D_n$ is used twice: once to construct the classifier and once for its evaluation. To control this optimism or reduce the bias, several strategies have been proposed, including cross-validation and bootstrapping. Their objective, in particular, is to estimate the risk without bias or with reduced bias by avoiding using the same data to construct the classifier and to study its performance. Some versions of these resampling methods are presented and also discussed according to the properties of the proposed risk estimators (the bias and the variance, $E[\hat{\mathcal{R}}]$ and $V[\hat{\mathcal{R}}]$, respectively, of an estimator $\hat{\mathcal{R}}$) in subsection 2.2.4.2.

Finally in subsection 2.2.4.3, we discuss the imbalance class problems for such classification purpose.

Note that this classification framework is used in the attribution method developed in Chapter 4 of this thesis.

### 2.2.4.1 The popular learning algorithms

We assume that the vector $X$ is composed of $p$ variables $X^1, X^2, \dots, X^p$ and recall that $Y$ is the label, i.e. the target variable.

#### Classification and Regression Trees (CART)

CART is a decision tree algorithm designed to create a tree-like structure by recursively partitioning the data. The root node of the tree contains all the data, each internal node corresponds to a split of the data into two
subsets and the terminal node, called leaves, is associated to a class label. The paths from root to leaf represent classification rules. An example of a classification tree is given in Figure 2.5 with comments in the caption.

<figure><img src="image_11.png" /><figcaption>Figure 2.5 – A tree showing survival of passengers on the Titanic ("sibsp" is the number of spouses or siblings aboard). The number in the leaves give the probability of survival and the percentage of observations. If the passenger is a male at most 9.5 years old with strictly fewer than 3 siblings, he will be classified as a survivor (Wikipedia contributors, 2023).</figcaption></figure>

The construction of the tree therefore consists in providing a sequence of nodes where each node is defined by a question, i.e. the joint choice of a variable $X^j$ among the $p$ variables and a threshold $d$ of the quantitative variable $X^j$:

$$ \{X^j \le d\} \cup \{X^j > d\}. $$

Splitting the data according to this decision rule means that the observations with a value of the $j$th variable smaller than $d$ go to the left child node, and all those with a value larger than $d$ go to the right child node (the two branches of the node). The objective is to make the best splitting, i.e. to search the best combination variable-threshold making the two resulting subsamples as pure/homogeneous as possible according to a criterion. The CART algorithm thus requires:

* a "purity" or "homogeneity" criterion. Common metrics include the Gini index, misclassification error, and entropy. In this study, we employ the Gini index given by $Gini = 1 - \sum_{k=1}^{K} p_k^2$ where $K$ is the number of classes and $p_k$ corresponds to the probability of an individual being classified to class $k$.

* a stopping rule since this process only halts when further splitting of a node is no longer feasible such as when each terminal node (leaf) contains only one observation.

* a leaf assignment. Lastly, each terminal node is assigned to a class label.

Choosing a stopping rule a priori is a challenging task. In practice, we construct a complete tree called the maximal tree and denoted $T_{max}$, that is to say a tree large enough to classify without error all the data in the training sample (leading inevitably to overfitting). Then we search for the best pruned subtree of $T_{max}$, i.e. from which we remove the branches that do not provide significant predictive power. An exhaustive search (considering all the subtrees) is not possible from an algorithmic point of view since their number is exponential. To get around this problem, Breiman et al. (1984) proposed to construct a sequence of nested
subtrees $T_1 \le T_2 \le \dots \le T_{\max}$ and to select the one that minimizes the following penalized misclassification error:

$$ \hat{\mathcal{R}}_n(T_i) + \alpha|T_i|, $$

where $|T_i|$ represents the size of the tree $T_i$ and $\hat{\mathcal{R}}_n(T_i)$ its empirical misclassification.

In conclusion, CART is a non-parametric method that is simple to understand, interpret, and visualize. The drawback is that it is particularly unstable, i.e. very sensitive to fluctuations in the sample (small variations in the data might result in a completely different decision tree).

### Random Forest (RF)

In order to overcome the issue of instability in CART, Breiman (1996) introduced the Random Forest (RF) algorithm. The idea of this algorithm is to consider multiple decision trees to improve the predictive power but at the cost of a loss of interpretability. According to the Strong Law of Large Numbers, the error rate of the ensemble of trees converges as the number of trees increases (Breiman, 2001). This convergence effectively mitigates the overfitting problem commonly associated with single decision trees, thereby enhancing both the model's accuracy and efficiency.

The underlying principle of RF can be succinctly explained by the Bagging (Bootstrap Aggregating) technique, which is a method for generating multiple versions of a predictor on bootstrap samples and using these to get an aggregated predictor (Breiman, 1996) (as illustrated in Figure 2.6). Further details on bootstrap sampling are provided in subsection 2.2.4.2, and the typical approach to result aggregation is by majority vote.

<figure><img src="image_12.png" /><figcaption>Figure 2.6 – Random forest illustration (TIBCO, 2023).</figcaption></figure>

In Random Forest, the decision rule (question) at each node is determined on the basis of a random subsample of variables $m$ out of $p$ variables. This additional layer of randomness serves to further reduce the correlation between the decision trees, reducing the variance of the aggregated predictions and thereby enhancing the bagging. To quantify the risk of the ensemble, Breiman (2001) introduced the concept of out-of-bag (OOB)
estimation, which is elaborated upon in subsection 2.2.4.2.

Random Forest offers several advantages, including enhanced predictive accuracy, and flexibility in ensemble composition, as it can incorporate various types of classifiers. However, it is essential to acknowledge that Random Forest can be computationally expensive, especially as the number of trees in the forest grows. Additionally, the aggregation of results from multiple trees can make the interpretation of the model more challenging compared to a single decision tree.

### k Nearest Neighbors (k-NN)

Another powerful non-parametric classification algorithm is the k-NN (Fix and Hodges, 1989; Cover and Hart, 1967). The intuitive and simple idea behind k-NN is to assume that very close observations probably belong to the same class. Figure 2.7 illustrates the k-NN procedure. Initially, we have data in classes A and B and need to classify an unknown observation, represented by the purple circle with a question mark, into one of these classes. The k-NN procedure involves the two following steps for classifying a new observation x:

1. Identify the *k* neighbors based on their distance to x (the distance metric can be the Euclidean distance, the Manhattan distance, or any other suitable measure). In the right plot of Figure 2.7, *k* = 3 neighbors are chosen.

2. Determine the predicted class of x by the majority class among the *k* identified neighbors. For this example, x is classified into class B.

<figure><img src="image_13.png" /><figcaption>Figure 2.7 – k-NN illustration (IBM, 2023).</figcaption></figure>

This algorithm is straightforward and does not require a separate training phase, making it relatively easy to implement. However, this method becomes computationally expensive when applied to large datasets. Computing the distance between observations becomes computationally complicated as the number of variables increases. Additionally, it tends to exhibit poorer performance on imbalanced datasets, where certain classes have substantially more observations than others. In such cases, k-NN may favor the majority class during prediction, leading to biased outcomes. The choice of *k* is of course crucial for the estimation quality: a small *k* may result in overfitting, while a large *k* may lead to underfitting. This number is classically chosen via cross-validation.
### Linear Discriminant Analysis (LDA)

Introduced in 1936 by Fisher (1936), Linear Discriminant Analysis (LDA) was initially designed for two-class problems. The objective of the LDA is to estimate, within a parametric framework, the Bayes classifier, $\psi^*$. More precisely, the rule of the Bayes classifier is to classify an observation $x$ in the most probable group, i.e. the group that maximizes the posterior probability of group membership given by

$$P(Y = k | X = x).$$

This conditional probability can be decomposed using the Bayes formula as follows:

$$P(Y = k | X = x) = \frac{\pi_k f_k(x)}{\sum_l \pi_l f_l(x)},$$

where the $\pi_k = P(Y = k)$ are the prior probabilities of group membership and the $f_k$ are the density distributions of $X | Y = k$.

The LDA assumes that the distribution of $X$ in each class is gaussian

$$X | Y = k \sim N(\mu_k, \Sigma),$$

where $\mu_k \in \mathbb{R}^p$ and $\Sigma$ is a definite positive matrix with size $p \times p$. It is straightforward to show that for two classes $k$ and $l$,

$$B(x) = \log \left( \frac{P(Y = k | X = x)}{P(Y = l | X = x)} \right) = \log \left( \frac{f_k(x) \pi_k}{f_l(x) \pi_l} \right) = b^T x + b_0,$$

with $b = \Sigma^{-1}(\mu_k - \mu_l)$ and $b_0 = -\frac{1}{2}(\mu_k^T \Sigma^{-1} \mu_k - \mu_l^T \Sigma^{-1} \mu_l) + \log(\pi_k / \pi_l)$. Between these two classes, an observation will be classified to class $k$ if $B(x) > 0$. The decision boundary between classes $k$ and $l$ is thus a linear function of $x$.

The prior probabilities $\pi_k$ and the distribution parameters $\mu_k$ and $\Sigma$ are estimated using the classical maximum likelihood method.

LDA offers multiple advantages such as computational efficiency, simplicity, and can work well even if the number of variables is large. Moreover, it can take into account multicollinearity (correlation between variables) in the data. The limitations of LDA mainly relate to the normality assumptions and the greater the deviation from this hypothesis, the less it will be possible to guarantee the performance of the chosen classifier. One can choose to consider that the covariance matrix is class-specific. In this case, the decision boundary is rather a quadratic function of $x$, called the QDA. The QDA is more flexible than the LDA but at the price of a high cost of the estimation.

### Parameters of the algorithms

All these algorithms involve parameters that the user must chosen (the penalty constant $\alpha$ in CART, the number of considered variables at each node $m$ for RF or the number of neighbors for kNN). Generally, these parameters are calibrated via cross-validation strategy applied on a predefined grid of values for the parameter. Then the
optimal value is the one that minimizes the estimated risk. In this thesis, we use the K = 10-fold cross-validation and the R package caret.

### 2.2.4.2 Resampling methods

As outlined in the introduction, the two most common used resampling methods to evaluate the performance of a classifier or model are cross-validation and bootstrap.

#### Cross validation

The general principle of cross-validation (CV) involves dividing the sample $D_n$ into two sub-samples: a training sample $D_n^L$ used to train a classifier and a testing sample $D_n^T$ (the remaining data, $D_n^T = (D_n^L)^c$ with size $n_T$) used to measure the performance of this classifier. Due to the independence between $D_n^L$ and $D_n^T$, we obtain a good assessment of the risk of $\hat{\psi}(D_n)$ (avoiding the over-optimism of the empirical risk) given by

$$ \hat{\mathcal{R}}^{\text{val}}(\hat{\psi}; D_n; L) = \frac{1}{n_T} \sum_{i \in T} \mathbb{1}_{\{\hat{\psi}(D_n^L, X_i) \neq Y_i\}} $$

Since various strategies can divide the same sample, there exists a large number of possible validation procedures. The common ones are:

* the "hold-out" validation or simple validation, which involves a single division of $D_n$. This procedure can be repeated K times, each time changing the validation sample. One limitation, however, is that this procedure does not guarantee consistent data dependencies in the test set across iterations.

* the "K-fold" cross-validation, which consists in dividing the sample $D_n$ into K blocks of observations of equal size, each block serving in turn as a validation set and the remaining K - 1 blocks composing the training sample. Note that the known leave-one-out cross-validation (LOOCV) is a particular case of the K-fold CV with K = n. In practice, K needs to be chosen (by the user) and this choice is crucial: a smaller K yields a higher bias of the risk estimator, but the variance is concurrently diminished. The well-known trade-off between bias and variance is detailed in Rodriguez et al. (2010).

The risk evaluated by CV is then the average of the K risks:

$$ \hat{\mathcal{R}}^{\text{CV}} = \frac{1}{K} \sum_{k=1}^{K} \hat{\mathcal{R}}^{\text{val}}(\hat{\psi}; D_n; L(k)) $$

where $L(k)$ is the kth set of observations belonging to the training set $D_n^L$, i.e. the training set of the kth repetition in the repeated hold-out validation or the training set including all blocks except the kth one (on which the constructed classifier is evaluated) in the K-fold CV.

As Arlot (2018) explains in his interesting document: "Intuitively, we can already say that K-fold cross-validation has the advantage of making a “balanced” use of the data: each observation is used exactly K - 1 times for training and once for learning. This is by no means guaranteed with repeated hold-out validation. On the other hand, one can wonder about the drawbacks of always using together (either for training or for validation) the
observations of the same block. The repeated hold-out approach, due to its random nature, makes it possible to avoid any biases induced by this link between observations".

### Bootstrap

In contrast to cross-validation procedures, which partition the original data to train and evaluate, the bootstrap method (Efron and Tibshirani, 1993) involves random sampling with replacement from the original dataset to create new datasets. These datasets, known as bootstrap samples of the same size, are used to construct classifiers, and then combined to evaluate the risk.

More precisely, the ordinary bootstrap procedure consists in repeating the following steps B times:

* create a new dataset, denoted $D_n^b$, from the original one, $D_n$, by randomly sampling $D_n$ with replacement,

* construct a classifier $\hat{\psi}^b = \hat{\psi}(D_n^b)$,

* calculate its risk on $D_n$: $\hat{\mathcal{R}}_n(\hat{\psi}^b; D_n)$ (where $\hat{\mathcal{R}}_n$ is given by (2.24)).

Combining all the B risks, the estimation of the risk by bootstrap is therefore:

$$ \hat{\mathcal{R}}^{\text{Boot}} = \frac{1}{B} \sum_{b=1}^{B} \hat{\mathcal{R}}_{n}(\hat{\psi}; D_{n}). $$

This estimator is generally biased due to optimism, but it does offer enhanced precision in the estimation of the risk.

To counter this optimistic bias, Breiman (2001) proposed the out-of-bag (OOB) bootstrap estimator. During sampling, certain data points are replicated within each "bag", leading to a portion of the original data being excluded and forming the "out-of-bag" set. The classifier derived from each sample is subsequently tested on its associated "out-of-bag" set. The risk estimator becomes:

$$ \hat{\mathcal{R}}^{\text{OOB-B}} = \frac{1}{n} \sum_{i=1}^{n} \frac{1}{B_i} \sum_{b \in K_i} \mathbb{1}_{\{\psi^b(X_i) \neq Y_i\}}. $$

where $K_i$ indicates the set of bootstrap samples that do not contain the observation i and $B_i$ is its size $B_i = |K_i|$.

### Summary

Both cross-validation and bootstrap methodologies are valuable for evaluating model performance, but they have different strengths and weaknesses. If the cross-validation is a relatively simple procedure to understand and implement, it can be computationally demanding especially for large datasets, and of course it is sensitive to the choice of the folds and their number. Bootstrap is a more flexible method but it is conceptually more difficult to understand and thus less used compared to CV. Note that for equal computational effort K = B, the bootstrap provides an estimator with a lower variance and the cross-validation with a lower bias.

In some scenarios, a combined approach can yield more robust error estimations. For instance, the leave-one-out bootstrap method, discussed by Efron and Tibshirani (1997) exhibits reduced variance compared to the
traditional leave-one-out cross-validation (LOOCV) approach.

Similarly, the concept of "bootstrap cross-validation", initially introduced by Fu et al. (2005), involves employing cross-validation within each bootstrap sample. This method offers multiple advantages, including an unbiased estimator of cross-validation for each individual sample and the ability to effectively handle imbalanced datasets through the utilization of the bootstrap strategy.

### 2.2.4.3 Imbalance class problem

One important problem in classification arises when the classes are imbalanced. The imbalance is often inherent to the problem at hand, as for example for medical diagnosis of a rare disease. In such situation, the learning algorithms will have difficulty to correctly identify the minority classes (the resulting classifier being strongly biased for the minority classes).

To address this problem, a typical approach involves a resampling technique that aims to rebalance the classes. Random oversampling of the minority class or undersampling of the majority class can be employed to modify the data and create a more balanced data set. While resampling can be a simple and straightforward solution that can improve classification performance, it can also introduce potential issues. Oversampling can lead to overfitting, as duplicate data is introduced, while undersampling can cause underfitting by removing representative data. In addition to resampling techniques, He and Garcia (2009) has reviewed other methods such as Cost-Sensitive methods, Kernel-based methods, and Active Learning methods. Among these methods, Cost-Sensitive methods prove to be particularly valuable as they take into account the cost of misclassifying instances in the minority class more than those in the majority class. This approach can be more effective than sampling methods, especially when the minority class carries greater significance. However, achieving optimal performance with Cost-Sensitive methods requires careful tuning of the cost parameters.
# Chapter 3

# Sensitivity of Change-Point Detection and Trend Estimates to GNSS IWV Time Series Properties

## 3.1 Summary

This chapter is presented in the form of a published article. The objective of the study is to investigate the sensitivity of the segmentation method to data properties and assess the consequential impact of these properties on trend estimations. We undertake this investigation by comparing pairs of datasets based on four critical factors: GNSS data processing methods (IGSrepro1 vs. CODE REPRO2015), the temporal extent of the time series (17 years vs. 25 years), the auxiliary data utilized in the Integrated Water Vapor (IWV) conversion, and the sources of reference data (ERAI and ERA5). Our analysis places particular emphasis on specific data properties, including mean values, noise levels, and periodic bias, as these properties have an impact on the segmentation results, particularly in terms of their influence on the number and positioning of change-points. This investigation allows us to both comprehend the segmentation results and identify important features of GNSS that can create inhomogeneities and impact trend estimates.

Our findings indicate significant impacts on the segmentation results when altering the GNSS processing and reference reanalysis. Only 45-49% of change-points are similar in these cases, compared to 71-81% similarity for the other two remaining factors (temporal extent and auxiliary data). Notable changes in GNSS processing include adjustments in the a priori ZHD correction, antenna/radome calibration model, and mapping function. Improvements in CODE processing lead to noise reduction and decreased periodic bias. Similarly, transitioning from ERAI to ERA5 as a reference reduces representativeness errors, resulting in noise and periodic bias reduction, making it easier to detect smaller change-points. The validation rate of the detected change-points with respect to metadata was found to consistently falls within the range of 30% to 35% for all datasets, i.e. the impact of data properties on that metric is small. This result suggests that the main instrumental changes that impact the GNSS IWV estimates were captured by the segmentation in both datasets.

The impact on trend estimations is investigated in two aspects. Firstly, we examine the influence of the four
factors previously explored in segmentation on IWV trend estimates. We find that changing the length of the time series has the strongest impact, both on the average and dispersion of the trend estimates across the network. The change in mean is believed to reflect the intensification of the water cycle in the 2010-2020 decade compared to the previous one (see the steeper global mean IWV trend in Figure 1.3). The reduced dispersion with the longer period is mainly a result of the smaller standard error of the estimates. Secondly, we consider the impact of inhomogeneity correction on trend estimates. We find that it is particularly pronounced when only change-points validated with respect to GNSS metadata (receiver, antenna, radome changes) are corrected. This impact is evident in global mean trends, dispersion, and RMS differences with respect to ERA5. Correction using all change-points also impacts the trend estimates but to a slightly lesser extent. This underscores the importance of addressing the attribution step, which aims at retaining only the change-points due to GNSS for the correction. The development of this step is the topic of Chapter 4. Finally, the dispersion associated with different homogenized trend estimates can be used as a measure of the uncertainty in the trend estimate at a single station. It amounts to 0.1 – 0.2 kg m⁻² decade⁻¹, or 0.5-1% decade⁻¹, which confirms the feasibility of detecting relevant global and regional climate trends with the GNSS IWV data.

## 3.2 Paper No. 1
Article

# Sensitivity of Change-Point Detection and Trend Estimates to GNSS IWV Time Series Properties

Khanh Ninh Nguyen<sup>1,2,*†</sup> ID, Annarosa Quarello<sup>1,2,†</sup>, Olivier Bock<sup>1,2,†</sup> ID and Emilie Lebarbier<sup>3,†</sup>

<figure><img src="image_14.png" /></figure>

check for
updates

**Citation:** Nguyen, K.N.; Quarello, A.; Bock, O.; Lebarbier, E. Sensitivity of Change-Point Detection and Trend Estimates to GNSS IWV Time Series Properties. *Atmosphere* **2021**, *12*, 1102. https://doi.org/10.3390/atmos12091102

Academic Editor: Peter Domonkos

Received: 30 July 2021
Accepted: 20 August 2021
Published: 26 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_15.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** This study investigates the sensitivity of the GNSSseg segmentation method to change in: GNSS data processing method, length of time series (17 to 25 years), auxiliary data used in the integrated water vapor (IWV) conversion, and reference time series used in the segmentation (ERA-Interim versus ERA5). Two GNSS data sets (IGS repro1 and CODE REPRO2015), representative of the first and second IGS reprocessing, were compared. Significant differences were found in the number and positions of detected change-points due to different a priori ZHD models, antenna/radome calibrations, and mapping functions. The more recent models used in the CODE solution have reduced noise and allow the segmentation to detect smaller offsets. Similarly, the more recent reanalysis ERA5 has reduced representativeness errors, improved quality compared to ERA-Interim, and achieves higher sensitivity of the segmentation. Only 45–50% of the detected change-points are similar between the two GNSS data sets or between the two reanalyses, compared to 70–80% when the length of the time series or the auxiliary data are changed. About 35% of the change-points are validated with respect to metadata. The uncertainty in the homogenized trends is estimated to be around 0.01–0.02 kg m⁻<fnref n="2" /> year⁻<fnref n="1" />.

**Keywords:** segmentation; homogenization; climate; GNSS; integrated water vapor; time series; trend; reanalysis

## 1. Introduction

Long records of observational data are essential to monitoring climate change and understanding the underlying climate processes [1,2]. However, long time series are often affected by inhomogeneities due to changes in instrumentation, in station location, in observation and processing methods, and/or in the measurement conditions around the station [3]. Inhomogeneities often take the form of abrupt changes, which are detrimental to estimating trends and multi-scale climate variability [4]. Various homogenization methods have been developed for the detection and correction of such change-points in the context of climate data analysis [1,2,5–9].

In this paper, we are interested in ground-based Global Navigation Satellite System (GNSS) integrated water vapor (IWV) measurements. GNSS measurements are qualified among the most accurate and continuous IWV measurements in all weather conditions but have only quite recently been considered for climate analysis [10–13]. Parracho et al. [14] was one among the first to analyze global IWV trends from more than 15 years of GNSS data and to confront them to the ECMWF reanalysis, ERA-Interim (ERAI), [15], and to the NASA/MERRA-2 reanalysis [16]. Significant differences were discovered in IWV trends

<footnote n="1">¹ Institut de Physique du Globe de Paris (IPGP), Centre National de la Recherche Scientifique (CNRS), Institut National de l'Information Géographique et Forestière (IGN), Université de Paris, 75005 Paris, France; annarosa.quarello@edu.unito.it (A.Q.); bock@ipgp.fr (Q.B.)</footnote>

<footnote n="2">² Ecole Nationale des Sciences Géographiques (ENSG), Institut National de L'information Géographique et Forestière (IGN), 77455 Marne-la-Vallée, France</footnote>

<footnote n="3">³ Laboratoire Modal'X, UPL, Université Paris Nanterre, 92000 Nanterre, France; emilie.lebarbier@parisnanterre.fr</footnote>

<footnote>* Correspondence: knguyen@ipgp.fr</footnote>

<footnote>† These authors contributed equally to this work.</footnote>
between the two reanalyses and between the reanalyses and the GNSS data. On the one hand, this study pointed to the importance of the atmospheric model, the assimilation system, but also the quality and quantity of assimilated observations in reanalyses. On the other hand, inhomogeneities were also suspected in the GNSS data at several sites. Developing a homogenized GNSS IWV time series is of prime importance to estimate regional and global IWV trends and variability but also to verify climate models and reanalyses. This study investigates in more detail the homogeneity of the GNSS IWV data set used by Parracho et al. [14], as well as a more recently reprocessed GNSS data set. It also updates the previous results from Parracho et al. [14] and Bock and Parracho [17] with the new ECMWF reanalysis, named ERA5 [18].

The main causes of inhomogeneities in GNSS IWV time series are:

* Equipment changes (antenna, radome, and receiver). Each antenna/radome pair has a particular impact on the measurements, which is taken into account at the processing level with a specific calibration model (see Section 2). However, model imperfections, multipath and on-site electromagnetic coupling with the antenna's environment, and equipment aging are responsible for small biases which can change over time. The quality of measurements also depends on the receivers. Modern receivers have more stable clocks, reduced cycle slips, and noise and are capable of observing satellites from new GNSS systems (GPS, GLONASS, etc.). Hence, changes in data quality/properties are expected, which can introduce offsets and possibly trends (e.g., when new satellites are introduced progressively). Changes in receiver settings, such as cutoff angle, are also known to produce abrupt changes in the mean IWV estimates [19].

* Changes in the environment near the receiver antenna can introduce multipath and obstructions that alter the measurements and cause inhomogeneities.

* Processing changes. The details of the data processing are known to impact the IWV estimates. The most important aspects and parameters are the tropospheric model (the mapping functions, the a priori hydrostatic model, the time-dependency), the antenna/radome calibration models, the elevation-dependent weighting, and the cutoff angle (see Section 2).

The first cause is well documented for International GNSS Service (IGS) stations and other scientific networks (ftp://igs.ign.fr/pub/igs/igscb/station/log/, accessed on 30 July 2021). Therefore, metadata can be used to check if change-points detected in the IWV time series can be explained by known equipment changes. The second cause is usually not well documented, but the analysis of the raw measurements and post-fit residuals can help to detect changes in the environment. The third cause is of a different nature as it depends on the analysis procedure and models, which are both the subject of active research in order to improve the accuracy and homogeneity of the GNSS products (see Section 2). However, not all biases and inhomogeneities can be corrected at the processing level, and further post-processing homogenization methods are needed.

Many different homogenization methods have been developed by climatologists. The heart of any homogenization method is the detection of change-points, the so-called segmentation method. Some segmentation methods use statistical tests [6,7], while others use a penalized likelihood approach [1,2,20]. The performance of both approaches are comparable, but, in general, the results depend on the data properties (nature of the background noise, presence of a periodic bias and/or a trend), the adopted model (parametric or non-parametric), and the search method (optimal or sub-optimal) [9,21].

Quarello [21] developed a segmentation method, called GNSSseg, especially devoted to detect changes in the mean of time series of IWV differences between GNSS and a reference and taking into account the presence of a periodic bias and a heterogeneous noise with a monthly variation. The method uses a penalized likelihood approach and is optimal in the sense that the estimation of the positions of the change-points is done using an efficient algorithm. The method proposes several penalty criteria, which aims to choose the number of change-points, with different sensitivities to the data properties (length of
the time series, noise distribution, etc.). The use of several criteria can help to mitigate their limitations but requires special post-processing to make the final decision, either automatic or manual. The post-processing may also include outlier detection, validation with metadata when available, and manual inspection. The automatic version of the GNSSseg algorithm was evaluated in a benchmark exercise and compared to other existing segmentation methods where it was found to be one of the most efficient in detecting change-points in synthetic time series mimicking the GNSS minus ERAI IWV differences at the moderate complexity [22].

The general objective of this paper is to evaluate the sensitivity of segmentation results with the GNSSseg method (recently improved in terms of computational time and so-called GNSSfast method) and the subsequent trend estimates to various qualitative and quantitative properties of both GNSS and reference data. The study considers the particular cases of two different GNSS data sets (IGS repro1 and CODE REPRO2015) combined with two different reanalysis data sets, ERA-Interim [15] and ERA5 [18], which serve as references to compute to IWV differences used in the segmentation. IGS repro1 and CODE REPRO2015 are representative of the 1st and 2nd generation of IGS reprocessing products, and, as such, they are expected to be of different quality. They also cover different time periods. ERAI and ERA5 are the 4th and 5th generation reanalyses produced by ECMWF [18] and are also of different quality and spatial resolution.

The paper is organized as follows. In Section 2, we describe the characteristics of the two GNSS IWV data sets and discuss which factors in the data processing control the accuracy of the daily IWV estimates and their homogeneity in the long term. We also present the global homogenization and the trend estimation methods. In Section 3.1, we study the impact of data properties on the segmentation results. The following questions are specifically investigated: (1) What is the impact of the different data processing between IGS repro1 and CODE REPRO2015 on the segmentation results? (2) what is the impact of the time length on the segmentation results? (3) What is the impact of the reference data source (ERAI versus ERA5)? (4) What is the impact of the auxiliary data source used in the ZTD to IWV conversion (ERAI versus ERA5)? The segmentation results will be compared on the basis of several statistics: the number of detections before and after outlier screening, the number of outliers, and the number of validations/attributions after screening using metadata and nearby stations. In Section 3.2, the corrected time series are used to estimate linear trends, and the same questions as above are investigated for the trends. Finally, Section 4 concludes the paper and discusses the future work.

## 2. Materials and Methods

### 2.1. GNSS IWV Data

Before describing the characteristics of the two GNSS IWV data sets used in this study, we deem it is necessary to discuss which factors in the data production control the accuracy of the daily IWV estimates and their homogeneity in the long term. The raw GNSS measurements consist of code and carrier phase signals transmitted by GNSS satellites (GPS, GLONASS, etc.), which are measured by ground-based stations composed of an antenna and a multi-channel receiver [23]. The measurements collected from a global tracking network are analyzed in routine by several International GNSS Service (IGS) analysis centers who compute accurate satellite orbits and clock parameters (https://www.igs.org/acc/, accessed on 30 July 2021). These, as well as other parameters (Earth Rotation Parameters (ERPs) and satellite phase biases), are necessary for subsequent users to process accurate station positions and tropospheric delays for their own applications (e.g., geodesy, surveying, weather forecasting, etc.). Both the IGS and user data processing procedures rely on the use of various models to account for geophysical effects (solid Earth tides, ocean tides), atmospheric propagation effects through the ionosphere and the neutral atmosphere, and instrumental biases and variations induced by the transmitter and receiver antennas and electronics [24]. The propagation effect in the neutral atmosphere, so-called tropospheric delay, is traditionally decomposed into its hydrostatic and wet components
in the zenith direction, which are mapped into the slant observation direction which their respective mapping functions, and a two-parameter horizontal gradient model [25]:

$$STD = ZHD \times m_h(\epsilon) + ZWD \times m_w(\epsilon) + (\vec{G} \cdot \vec{e}) \times m_G(\epsilon), \quad (1)$$

where ZHD is the Zenith Hydrostatic Delay (ZHD), mainly due to the contribution from dry air, ZWD is the Zenith Wet Delay (ZWD), due to water vapor molecules, $\vec{G}$ is the horizontal gradient vector, $\vec{e}$ is the unit vector pointing into the direction of the observed satellite, $m_h$, $m_w$, and $m_G$ are the respective mapping functions, and $\epsilon$ is the elevation angle.

During the GNSS data processing, ZHD is corrected a priori, while ZWD and G are estimated. The ZWD estimates can include some residual bias from incorrect a priori ZHD correction or a deficiency of the hydrostatic and wet mapping functions. The GNSS software traditionally provide the total Zenith Tropospheric Delay (ZTD), which is the sum of the a priori ZHD, ZHD<sub>ap</sub>, and estimated ZWD, ZWD<sub>est</sub>:

$$ZTD = ZHD_{ap} + ZWD_{est}. \quad (2)$$

The main parameter of interest for climate studies is the Integrated Water Vapor (IWV), i.e., the integral of the water vapor molecules at the zenith. The conversion of ZTD into IWV operates, thus, in two steps [25]:

$$ZWD = ZTD - ZHD \quad \text{and} \quad \text{IWV} = \kappa(T_m) \times ZWD. \quad (3)$$

The accuracy of the IWV estimates is, thus, basically determined by the accuracy of the ZTD parameters derived from the GNSS data processing and the quality of the ZWD to IWV conversion procedure. For the purpose of retrieving IWV, a more accurate estimate of ZHD is required than the one used a priori for the processing. It can be derived from the surface air pressure data, $P_s$, available at the GNSS station [26]:

$$ZHD = 10^{-6} k_1 R_d \frac{P_s}{g_m}, \quad (4)$$

where $k_1$ is the dry air refractivity coefficient, $R_d$ is the dry air specific gas constant, $P_s$ is the surface air pressure, and $g_m$ is the mean acceleration due to gravity. The ZWD to IWV conversion factor $\kappa$ is defined as a function of the weighted mean temperature $T_m$ [10]:

$$\kappa(T_m) = \frac{10^6}{R_v(k'_2 + \frac{k_3}{T_m})},$$

where $k'_2$ and $k_3$ are refractivity coefficients for the water molecule, and $R_v$ is the specific gas constant for water vapor. The refractivity coefficients were taken from Bock [27]. The weighted mean temperature is defined by Bevis et al. [10]:

$$T_m = \frac{\int \rho_v(z) dz}{\int \frac{\rho_v(z)}{T(z)} dz}, \quad (5)$$

where $\rho_v(z)$ and $T(z)$ are the specific mass of water vapor and the air temperature, respectively, at height z above the surface. The integral is from the surface to the top of the atmosphere. It can be computed from a vertical profile of $\rho_v(z)$ and $T(z)$ given by a radiosonde climatology or an atmospheric model.

In this study, the auxiliary data, $P_s$ and $T_m$, required for the conversion of the ZTD estimates into IWV are computed from a global atmospheric reanalysis. Using reanalysis data has several advantages: the data are available at any position and time on the globe, the pressure and temperature are well constrained by observations, making the reanalysis data the best estimate of the global atmospheric state at any position and time, and the assimilation system uses an efficient screening and bias correction procedure to reject
suspect observational data and adjust bias changes associated to observational system
changes (e.g., between older and new satellites). In this study, we will consider two
different reanalyses: ERA-Interim [15] and ERA5 [18]. Details on the reanalyses are given
in the subsection below. The *P*<sub>s</sub> and *T*<sub>m</sub> data are computed from 6-hourly pressure level
data, as described in Reference [14] and the subsequent IWV estimates are aggregated into
daily estimates, as described in Reference [17]. The interest of comparing two different
reanalysis is that ERA5 is of superior quality due to the assimilation of more observations
and of higher spatial resolution, i.e., providing more accurate pressure and temperature
estimates in regions of steep topography. The accuracy of reanalysis estimates of *P*<sub>s</sub> and *T*<sub>m</sub>,
as compared to other data sources, is further discussed in Bock [27] and Bock et al. [28].

The main factors conditioning the accuracy of the ZTD estimates at the GNSS data
processing level are, by decreasing order, the hydrostatic and wet mapping functions,
the a priori ZHD correction data, and the antenna phase center variation (PCV) models.
Any bias in these data and models would map directly into the ZWD estimates. The
mapping functions are the most important because they determine how the signals from the
satellites at various elevations are mapped into the zenith direction (from a mathematical
point of view, they represent the partial derivatives, or regressors, of the ZHD and ZWD
parameters). Because the hydrostatic delay is corrected a priori, any bias in the a priori
ZHD, or in the hydrostatic mapping function, will also map into the estimated ZWD [29].
Receiver antenna phase center variations also highly depend on the elevation angle and, to
a lesser extent, on the azimuth angle. Before 2005, relative calibration models were used
in the IGS network, which were progressively replaced with absolute calibrations [30].
The convention at IGS is to use a type-mean calibration, i.e., a mean calibration model
determined from several calibrated antennas of the same type (producer and product)
when available, and a specific variant for each antenna/radome combination. When no
specific antenna/radome combination exists, the rule is to “adopt” the antenna calibration
without radome. In addition, when the antenna calibration does not exist, the calibration is
“copied” from a similar antenna (based on electronic and mechanic properties). The official
antenna/radome calibrations are distributed by IGS in the form of ANTEX (Antenna
Exchange Format) files, e.g., igs05.atx and igs08.atx, at the times when the IGS and CODE
data sets used in this study were produced. Absolute calibrations from “Robot” and
“Chamber” are the most accurate, while relative calibrations (type “field”) and relative
converted to absolute (type “converted”) are the less accurate [31]. The impact of satellite
and receiver antenna offsets (PCO) and PCVs on geodetic parameters have been mainly
investigated for positioning purposes. The impact on ZTD estimates has not been much
studied yet, although it is known that it would be tightly correlated with the vertical
position component. One of the goals of this paper is, thus, to examine how the change
from igs05.atx to igs08.atx impacts the accuracy and homogeneity of the GNSS IWV series
and to which extent these differences are detected by our segmentation method. In addition
to the aforementioned factors, there are other errors sources that can impact the accuracy of
the ZTD estimates, such as multipath and ambiguity fixing errors, as well as satellite orbits
and clock errors, and unmodeled and mis-modeled station displacements at sub-daily time
scales (e.g., tides). However, these are minor errors sources for the purpose of our study
here. For further details, see the study of Ning et al. [32]. Of main concern here are the
sources of bias and the mechanisms through which these biases can change with time, i.e.,
translate into inhomogeneities.

In this study, we consider two different GNSS data sets, which are representative of
two different generations of reprocessing products delivered by IGS (see Table 1). The first
one, referred to as IGS repro1, was produced by JPL/NASA in 2010–2011 in precise point
positioning (PPP) mode with GIPSY OASIS II software [33] as a special release of ZTD
estimates. This data set used the reprocessed IGS orbits, clocks, and ERPs produced by
JPL/NASA in the framework of the 1st reprocessing campaign organized by IGS. The repro-
cessed satellite products were generated for the period 1 January 1995 to 31 December 2007,
but JPL completed the series until mid-2011 in a consistent way. In this study, we use the
ZTD estimates until 31 December 2010 to have an integer number of years. According to the discussion above, the prominent features of the processing procedure are:

* Standard Temperature and Pressure (STP) model used for a priori ZHD correction [34],

* Global mapping function (GMF) for the hydrostatic and wet delays [35],

* IGS05 reference frame and igs05.atx absolute antenna PCO/PCV models

* 7° elevation cutoff angle, and

* uniform observation weighting.

<table><caption>**Table 1.** Processing strategies used by JPL/NASA to produce the IGS first reprocessed tropospheric data set (IGS repro1) and by the Center for Orbit Determination in Europe (CODE) to produce the CODE REPRO2015 reprocessed data set.</caption><thead><tr><th></th><th>IGS Repro1</th><th>CODE REPRO2015</th></tr></thead><tbody><tr><td>Software</td><td>GIPSY OASIS II</td><td>Bernese GNSS software v5.3</td></tr><tr><td>Strategy</td><td>PPP solution</td><td>Double-difference solution of a global network</td></tr><tr><td>Orbits, clocks, ERPs</td><td>IGS repro1 (1995.0–2008.0) + IGS final (2008.0–2011.0)</td><td>CODE repro2 (1994.0–2015.0) + CODE final (2015.0–2019.0)</td></tr><tr><td>Reference frame</td><td>IGS05</td><td>IGb08</td></tr><tr><td>Antenna calibration</td><td>igs05.atx</td><td>igs08_1852.atx until 28 January 2017, igs14.atx from 29 January 2017</td></tr><tr><td>Window length</td><td>24 h</td><td>72 h</td></tr><tr><td>Elevation cutoff angle</td><td>7°</td><td>3°</td></tr><tr><td>Observations</td><td>GPS</td><td>GPS (1994.0–2002.0), GPS+GLONASS (2002.0–2019.0)</td></tr><tr><td>Observation sampling</td><td>5 min</td><td>3 min</td></tr><tr><td>Observation weighting</td><td>uniform</td><td>σ² = 36 × 10⁻⁶ / cos²(Z) where Z = zenith angle</td></tr><tr><td>Tropospheric model</td><td>ZHD and ZWD a priori:<br>ZHD = 1.013 × 2.27 × exp(-0.116*ht), ZWD = 0.1 m.<br>GMF mapping functions (hydrostatic and wet).<br>Random Walk model for ZWD and gradient parameters with constraints:<br>3 mm h⁻¹/² (ZWD) and 0.3 mm h⁻¹/² (gradients).<br>ZWD and gradient sampling: 5 min</td><td>ZHD and ZWD a priori:<br>6-hourly ECMWF analysis (provided by TUV).<br>VMF1 mapping functions (hydrostatic and wet).<br>Piece-wise linear model for ZWD with constraints:<br>5 m absolute and 5 m relative.<br>Sampling: 2 h (ZWD), 24 h (gradients).</td></tr><tr><td>Tropo files</td><td>ZTD and gradient estimates provided in SINEX files (0000, 0005, ... 2345 UTC)</td><td>ZTD and gradient estimates provided in SINEX files (resampled to 01, 03, ... 23 UTC)</td></tr><tr><td>Coordinate estimates</td><td>Estimated once per 24 h</td><td>Estimated once per 24 h</td></tr><tr><td>Ambiguity resolution</td><td>Float</td><td>Fixed</td></tr></tbody></table>

It should be noted that, due to a flaw in the data file handling, the ZTD files for a number of stations are not available on the FTP repository, and older files were used instead, for which a major difference was the use of the older NMF hydrostatic and wet mapping functions [36].

The second GNSS data set, referred to as CODE REPRO2015, was processed by the Center for Orbit Determination in Europe (CODE) in 2015 [37] using the Bernese GNSS software [38] in network mode (double-differenced observations). This data set used the reprocessed IGS orbits, clocks, and ERPs produced by CODE in the framework of the 2nd reprocessing campaign organized by IGS. The reprocessed products cover the period from 1 January 1994 to 31 December 2014. For the purpose of this study, we completed the series with the operational CODE products until the end of 2018 [39]. However, the latter underwent a switch in the reference frame to IGS14 and in the antenna cnetwork alibration (double-differenced observations).model to igs14.atx, on 29 January 2017. This change was, thus, included in the metadata for the CODE solution.
The prominent features of the CODE processing procedure are:

* ECWMF grid estimates of a priori ZHD and ZWD distributed by Technical University of Vienna [40],

* Vienna mapping function (VMF1) for the hydrostatic and wet delays [40],

* IGb08 reference frame and igs08_1852.atx absolute antenna PCO/PCV models,

* 3° elevation cutoff angle,

* $1/\cos(zenith)^2$ observation weighting, and

* GPS observations (1994–2001) and GPS + GLONASS (2002–2014).

According to the discussion above, we see that all the main factors conditioning the accuracy of the ZTD estimates are different between the two data sets. Not discussed above are the elevation cutoff angle and observation weighting, which contribute significantly to the sensitivity of ZTD estimates to the various elevation-dependent error sources (mapping functions, antenna PCO/PCV models, and multipath). The purpose of using a lower elevation cutoff angle in CODE is to include more observations, i.e., increase the precision of the estimated parameters. However, multipath is generally higher at low elevations. To mitigate it, the lower elevation observations are down-weighted. The JPL/NASA processing strategy was different as they used a 7° cutoff angle and no down-weighting. Possibly, this strategy might be more sensitive to multipath and anomalous propagation effects at low elevations.

The ZTD estimates from both data sets were screened for outliers following the methodology described in Bock [41] and Stepniak et al. [42], and converted to IWV using either ERA-Interim or ERA5 reanalysis. The 6-hourly IWV data were compared to the reanalysis IWV data and further screened for outliers (for each station, the IWV differences exceeding the median ± five standard deviations were removed). Afterward, the IWV values from GNSS and reanalysis, and the IWV differences between GNSS and reanalysis, were aggregated into daily and monthly estimates and made publicly available on the AERIS data center [43,44].

Figure 1 shows the location of the GNSS station available from the two data sets. In this study, we selected 81 common stations, for which the time series in both data sets covered a period of at least 15 years.

<figure><img src="image_16.png" /><figcaption>**Figure 1.** Map of the GNSS stations available from the two reprocessed data sets: IGS repro1 (empty circles), CODE REPRO2015 (small dots), and the 81 common stations (full circles) used in this study.</figcaption></figure>

## 2.2. Reference IWV Data

Our homogenization method operates on IWV differences between a GNSS series and a reference series. Because the IGS network is quite sparse, we cannot use a nearby station as is commonly done by climatologists (as in Venema et al. [9]). Instead, for every
GNSS station, a series of IWV from each of the two reanalyses is derived, and daily IWV differences are formed, as explained above. In earlier studies, we found that ERA-Interim and GNSS IWV had significant representativeness differences in Antarctica and in regions of steep topography (Andes, Himalayas, etc.) or near the oceans [17]. In this study, we will investigate the impact of representativeness errors on the segmentation results by comparing the results from the two reanalyses. The spatial resolution of the reanalyses is 0.75° × 0.75° for ERA-Interim and 0.25° × 0.25° for ERA5. Reduced representativeness errors are, thus, expected from ERA5 data. Moreover, the IWV values computed from ERA5 are also expected to be of higher quality since this reanalysis used a more recent model and assimilation system, and assimilated a much larger number of observations, especially in recent years [18].

## 2.3. Homogenization Method

Figure 2 shows the data flow chart starting with the GNSS ZTD data and ending with the corrected IWV series. The first two steps (Conversion and Comparison) are described in the previous sub-section.

<figure><img src="image_17.png" /><figcaption>**Figure 2.** Flowchart of the general homogenization procedure.</figcaption></figure>

The third step is the segmentation, i.e., the detection of change-points in the mean of the IWV difference time series. Here, we use the fast version of the GNSSfast R package published by Quarello et al. [45]. This version is available on https://github.com/arq16/GNSSfast.git, accessed on 30 July 2021. The segmentation method estimates K, the
number of segments, or the number of change-points K-1, the K means of the segments, and the K - 1 positions of the change-points, as well as a periodic bias function f and the monthly variance of the noise (also with annual periodicity), in a unified model for each IWV difference time series. It has been shown in a previous paper that modeling a periodic bias function and a heterogeneous noise helps to accommodate for the representativeness differences between the GNSS IWV data and the reanalysis IWV data (see Quarello [21]). The inference procedure is based on a penalized likelihood approach. Several penalty criteria are considered, but here we will present only the results for the so-called Birgé and Massart's penalty based on the 'dimension jump' proposed by Birgé and Massart [46] and calibrated in Lebarbier [47]. In a previous study, the GNSSseg method was compared to other existing segmentation methods on the basis of a benchmark data set and was found to be one of the most efficient [22].

The noise in the IWV differences is generally not perfectly Gaussian and strong peaks can lead to spurious change-points. These false detections are checked and screened out by a special post-processing procedure, which is symbolized in Figure 2 by the Screening step. In this step, clusters of 2 or more nearby change-points are checked for a significant change in the mean before and after using a weighted t-test. In this study, we consider as a cluster all the change-points closer than 80 days (the value of 80 days was found as optimal based on a mixture model analysis). The change-points within a cluster are flagged as "outliers" in the following. If the change in the mean is not significant (at the level of 5%), all the change-points in the cluster are removed. If it is significant, the group of change-points is replaced with one change-point, and its position is taken as the mean of the positions.

With a relative segmentation approach, such as used here, climatologists often assume that the reference is homogeneous [9]. If the reference is not homogeneous, multiple comparisons are necessary to determine in a statistical way to which of the series each detected change-points belong to. This task is accomplished during the Attribution/Validation step. The attribution of a change-point to GNSS or reference is decided based on the comparison with nearby GNSS series, while the validation is referring to the comparison of the detected change-points with the GNSS metadata. For the attribution step, the idea is to use all the available stations in the GNSS data, as long as the time series covers a period of time encompassing the tested change-point. We tested the procedure with the CODE REPRO2015 data set, including all 434 stations, with a maximum distance of 200 km and a time window of ±six months around the change-points. Based on the available data, about 30% of the change-points from the selected 81 stations could be checked. This number is too small to apply this test in a systematic way to all the data sets used in this study. So, we decided not to include it in the general discussion but only for a few cases studies. Instead, we will only use the validation step with respect to the GNSS metadata with a window of ±62 days, as in Van Malderen et al. [22].

The last step is the correction of the GNSS time series for the changes in the mean based on the attributed/validated change-points. Here, we follow the approach of Van Malderen et al. [22] which consists in subtracting from the GNSS IWV series a piece-wise constant function constructed from the change in means of the GNSS-reference series, where the last segment is taken as a reference. The performance of the correction depends on all the previous steps.

## 2.4. Trend Estimation Method

Following Weatherhead et al. [48], we use a linear trend model:

$$Y_t = \mu + \omega X_t + S_t + N_t, \qquad (6)$$

where $Y_t$ is the IWV time series, $X_t$ the linear trend function, $S_t$ a seasonal component which will be represented by a fourth order Fourier Series, $S_t = \sum_{i=1}^{4} a_i \cos(2\pi i t/T) + b_i \sin(2\pi i t/T)$, t is time in days since some reference date, T = 365.25 days, and $N_t$ the noise which is assumed to be autoregressive of the order 1 (AR(1)), that is, $N_t = \phi N_{t-1} + \epsilon_t$, where the $\epsilon_t$ are independent random variables with zero mean and variance $\sigma_e^2$. The AR(1)