Analysis of global trends of total column water vapour from multiple years of OMI observations

In this study, we investigate trends in total column water vapour (TCWV) retrieved from measurements of the Ozone Monitoring Instrument (OMI) for the time range between January 2005 to December 2020. The trend analysis reveals, on global average, an annual increase in the TCWV amount of approximately +0.054 kg m⁻² yr⁻¹ or +0.21 % yr⁻¹. After the application of a Z test (to the significance level of 5 %) and a false discovery rate (FDR) test to the results of the trend analysis, mainly positive trends remain, in particular over the northern subtropics in the eastern Pacific.

Combining the relative TCWV trends with trends in air temperature, we also analyse trends in relative humidity (RH) on the local scale. This analysis reveals that the assumption of temporally invariant RH is not always fulfilled, as we obtain increasing and decreasing RH trends over large areas of the ocean and land surface and also observe that these trends are not limited to arid and humid regions, respectively. For instance, we find decreasing RH trends over the (humid) tropical Pacific Ocean in the region of the Intertropical Convergence Zone. Interestingly, these decreasing RH trends in the tropical Pacific Ocean coincide well with decreasing trends in precipitation.

Moreover, by combining the trends of TCWV, surface temperature, and precipitation, we derive trends for the global water vapour turnover time (TUT) of approximately +0.02 d yr⁻¹. Also, we obtain a TUT rate of change of around 8.4 % K⁻¹, which is 2 to 3 times higher than the values obtained in previous studies.

1 Introduction

Water vapour is the most abundant greenhouse gas in the Earth's atmosphere and is involved in several atmospheric processes across all atmospheric scales, starting from phenomena like cloud droplet growth on the microscale, thunderstorms on the mesoscale, and hurricanes on the synoptic scale and finally to the climate or global scale by influencing the Earth's energy balance via the greenhouse effect and cloud, lapse rate, and water vapour feedback mechanisms (Kiehl and Trenberth, 1997; Randall et al., 2007). According to the Clausius–Clapeyron (CC) equation, changes in saturated water vapour are closely linked to changes in air temperature, as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}E}{E}}={\displaystyle \frac{L_{\mathrm{v}}\left(T\right)}{R_{\mathrm{v}}}}{\displaystyle \frac{\mathrm{d}T}{T^{\mathrm{2}}}},\end{array}$$

with saturation water vapour pressure E, latent heat of vaporisation L_v, the specific heat capacity of water vapour R_v, and the air temperature T. For typical atmospheric conditions, the CC equation yields that, for a temperature increase of 1 K, it can be expected that the water vapour concentration increases by approximately 6 %–7 % if the relative humidity remains unchanged (Held and Soden, 2000). Thus, given its key role in many atmospheric processes, and considering the global warming of the atmosphere and ocean within the last few decades, accurate monitoring of changes in the global water vapour distribution is essential not only for a better understanding of the Earth's hydrological cycle but also of the climate system in general.

Several quantities exist to characterise the content of water vapour in the atmosphere. To determine the distribution of these quantities on the global scale, satellite missions offer great opportunities. Depending on the spectral range, satellite instruments can provide different information. For example, in the radio and thermal infrared spectral range it is possible to retrieve information of the vertical profile of the water vapour concentration (e.g. Kursinski et al., 1997; Susskind et al., 2003). Another important quantity is the water vapour content integrated over the complete atmospheric column, also known as integrated water vapour or total column water vapour (TCWV). In addition to the spectral ranges already mentioned, this quantity can be retrieved in the microwave (Rosenkranz, 2001; Wentz, 2015), in the shortwave and near-infrared (Bennartz and Fischer, 2001; Gao and Kaufman, 2003), and in the visible spectral range (e.g. Noël et al., 1999; Lang et al., 2003; Wagner et al., 2003; Grossi et al., 2015; Borger et al., 2020).

Based on these satellite observations, several studies in the past have investigated trends or changes in the global water vapour distribution (e.g. Trenberth et al., 2005; Wagner et al., 2006; Mieruch et al., 2008; Wang et al., 2016) and found rates of change that correspond to the CC response (e.g. Trenberth et al., 2005). Trenberth et al. (2005) analysed trends for the time period of 1988 to 2003 from a TCWV data set of merged microwave satellite sensors and found generally positive trends that are consistent with assumption of fairly constant relative humidity. Mieruch et al. (2008) combined TCWV measurements from GOME and SCIAMACHY in the visible red spectral range and also determined positive TCWV trends for the time period January 1996 to December 2003. More recently, Wang et al. (2016) investigated TCWV trends for the time period from 1995 to 2011 for a TCWV data set combining measurements from radiosondes, GPS radio occultation, and microwave satellite instruments. They found positive but slightly weaker TCWV trends which they attributed to the slowdown in the global warming rate since 2000 that had terminated in 2014.

Nevertheless, a major limitation of the assumption of a CC response is the assumption of temporally invariant relative humidity. Typically, it is assumed that the relative humidity close to the surface (especially over the ocean) remains constant, which was also confirmed by Dai (2006). Over land surfaces, however, this assumption is not always given (Simmons et al., 2010; Fasullo, 2012). For instance, Dunn et al. (2017) showed, with their observational data, that there has been first a constant and then a clear decrease in near-surface relative humidity over land masses since 2000.

In this study, we continue the analysis of the trends in TCWV. For this purpose, we are using an observational TCWV data set (Borger et al., 2021a) based on measurements of the Ozone Monitoring Instrument (OMI; Levelt et al., 2006, 2018) in the visible blue spectral range. In doing so, we investigate not only how strong the trends in water vapour are on the local scale but also to what extent the assumption of constant relative humidity is fulfilled there. Moreover, we also investigate how sensitive the global atmospheric water cycle (more specifically, the water vapour residence time) responds to changes in surface air temperature.

For this purpose, the paper is structured as follows: in Sect. 2, we briefly introduce the OMI TCWV data set and describe the scheme for the trend analysis in detail. Then, in Sect. 3, we present the trend results from the OMI TCWV data set and put these results in context to the trend results from other data sets. In Sect. 4, we analyse local trends in relative humidity derived from the OMI TCWV trends and investigate how these are related to changes in precipitation in Sect. 5. Moreover, in Sect. 6, we analyse the responses of the water vapour residence time to global warming. Finally, in Sect. 7, we will briefly summarise our results and draw conclusions.

2  Data set and methodology

2.1  MPIC OMI TCWV data set

For our study, we use the monthly mean MPIC OMI TCWV data set from Borger et al. (2021a, b). The data set is based on measurements of the Ozone Monitoring Instrument (OMI; Levelt et al., 2006, 2018), which are analysed by means of differential optical absorption spectroscopy (DOAS; Platt and Stutz, 2008) in the visible blue spectral range, using the TROPOMI (TROPOspheric Monitoring Instrument) TCWV retrieval of Borger et al. (2020). First, a spectral analysis is performed in a fit window of 430–450 nm, taking into account the specific instrumental properties of OMI (more details in Borger et al., 2021a). Then, these fit results are converted to TCWV via an iterative algorithm that finds the optimal water vapour profile shape.

The data set covers the time period from January 2005 to December 2020 and provides the TCWV values on a spatial resolution of ${\mathrm{1}}^{\circ }\times {\mathrm{1}}^{\circ }$. In an extensive validation study, Borger et al. (2021a) showed that the data set is in good overall agreement to other reference data sets such as RSS SSM/I (Mears et al., 2015; Wentz, 2015) or ERA5 (Hersbach et al., 2020), especially over ocean surface. Moreover, Borger et al. (2021a) demonstrated, in a temporal stability analysis, that their data set is consistent with the temporal changes in the reference data sets, and that it shows no significant deviation trends (i.e. relative deviation trends smaller than 1 % per decade), which is particularly important for climate studies.

The major advantages of this TCWV data set in comparison to others are that, on the one hand, the data set provides a consistent time series since it is based on measurements from only one satellite instrument. Thus, inter-instrumental offsets do not have to be corrected when merging the data time series of the different instruments. On the other hand, in contrast to other spectral ranges, TCWV retrievals in the visible blue spectral range have a similar sensitivity over ocean and land surfaces and thus allow for consistent global analyses.

2.2  Trend analysis

For the trend analysis, we follow the approaches of Weatherhead et al. (1998), Mieruch et al. (2008), and Schröder et al. (2016), in which the fit function is given as follows:

$$\begin{array}{}\text{(2)}& Y_t=m+b\cdot X_t+S_t+{\mathrm{\Theta }}_t+N_t={\mathbf{M}}_{t}x+N_t,\end{array}$$

with the intercept m, the slope or trend b, respectively, the increasing time index X_t, the seasonal components S_t, and a component accounting for the influence of geophysical teleconnections (e.g. the El Niño–Southern Oscillation, ENSO), Θ_t, which can all be summarised in a matrix M_t. The term N_t stands for the fit residuals with respect to the measurement time series.

The seasonal components are modelled as a sum of sine and cosine functions with up to four frequencies, as follows:

$$\begin{array}{}\text{(3)}& S_t=\sum _{i=\mathrm{1}}^{\mathrm{4}}\left[c_i\sin (i\cdot \mathit{\omega }{X}_t)+d_i\cos (i\cdot \mathit{\omega }{X}_t)\right],\end{array}$$

with $\mathit{\omega }=\frac{\mathrm{2}\mathit{\pi }}{\mathrm{12}}$.

To account for the influence of teleconnections, we include several teleconnection indices Ω_i in the trend analysis. For the case of ENSO, we include the NOAA Oceanic Niño Index (ONI), which, according to Wagner et al. (2021), has the strongest impact on the TCWV time series distribution. Moreover, we follow the recommendations from Trenberth and Stepaniak (2001) and include a second ENSO index. In our case, we apply the Trans-Niño Index (TNI; Trenberth and Stepaniak, 2001). Furthermore, we investigated the influence of several other teleconnection indices and found that the Pacific Meridional Mode (PMM) sea surface temperature index (Chiang and Vimont, 2004) has a particularly strong influence on the autocorrelation of the noise in the Pacific Ocean. Typically, trends are already removed from teleconnection indices. However, since the time series of the indices cover several decades, the detrending is optimised for this large time period. Accordingly, we have detrended the indices again for our chosen time period (2005–2020). Apart from the three detrended index time series themselves, their derivatives are also considered within the trend analysis, as follows:

$$\begin{array}{}\text{(4)}& {\mathrm{\Theta }}_t=\sum _{i=\mathrm{1}}^{\mathrm{3}}{\mathit{\theta }}_{\mathrm{1},i}\cdot {\mathrm{\Omega }}_i+{\mathit{\theta }}_{\mathrm{2},i}\cdot {\displaystyle \frac{\partial {\mathrm{\Omega }}_i}{\partial t}}.\end{array}$$

For the fit residuals N_t, we assume that they follow a first-order autoregressive process AR(1), which can be described as follows:

$$\begin{array}{}\text{(5)}& N_t=\mathit{\varphi }{N}_{t-\mathrm{1}}+{\mathit{\epsilon }}_t,\end{array}$$

with the autocorrelation ϕ. In classical statistical methods it is often assumed that data are independent. However, this is not always the case in environmental data, in particular for time series analysis, in which data are likely temporally autocorrelated. Thus, not accounting for autocorrelation can give misleading results when these classical statistical test methods are applied to strongly persistent time series (von Storch, 1999; Wilks, 2011). For instance, Weatherhead et al. (1998) showed that, in the presence of temporal autocorrelation, the uncertainty of a linear trend is linked to the level of autocorrelation as follows:

$$\begin{array}{}\text{(6)}& {\mathit{\sigma }}_{\mathrm{trend}}^{\mathrm{2}}\propto {\mathit{\sigma }}_N^{\mathrm{2}}\cdot {\displaystyle \frac{\mathrm{1}+\mathit{\varphi }}{\mathrm{1}-\mathit{\varphi }}}\propto {\displaystyle \frac{{\mathit{\sigma }}_{\mathit{\epsilon }}^{\mathrm{2}}}{\mathrm{1}-{\mathit{\varphi }}^{\mathrm{2}}}}\cdot {\displaystyle \frac{\mathrm{1}+\mathit{\varphi }}{\mathrm{1}-\mathit{\varphi }}},\end{array}$$

with the fit error ${\mathit{\sigma }}_N^{\mathrm{2}}$ influenced by the autocorrelation and the “true” fit error ${\mathit{\sigma }}_{\mathit{\epsilon }}^{\mathrm{2}}$. Consequently, positive (negative) autocorrelation can lead to an underestimation (overestimation) of the uncertainty of the trend which, in turn, can cause misleading results when classical statistical test methods (e.g. Z test) are used to classify if a trend is significant or not. Moreover, as the fit is not statistically efficient (i.e. it does not have the minimal variance), the fit results can also deviate from the “truth” (see also Appendix A).

Hence, to account for the effect of autocorrelation, we use the Prais–Winsten transformation (Prais and Winsten, 1954) and proceed as follows. First, to calculate the autocorrelation ϕ of the residuals, we perform a linear least squares fit of Eq. (2) to the time series of the TCWV data set as the first guess for each grid cell which yields the time series of N_t. Then, we estimate the autocorrelation function using the Gaussian-kernel-based cross-correlation function algorithm, as described in Rehfeld et al. (2011), via the NEST package (http://tocsy.pik-potsdam.de/nest.php, last access: 7 June 2022). The advantage of this algorithm is that it takes into account the complete data of an irregular spaced time series. From the autocorrelation function, the lag 1 autocorrelation ϕ can then be derived by simple linear algebra.

Figure 1 Global distribution of the lag 1 autocorrelation coefficients of the fit residuals (or fit noise) of the trend analysis for the MPIC OMI TCWV data set.

Figure 1 illustrates the global distribution of the lag 1 autocorrelation coefficients of the fit residuals from the trend analysis of the OMI TCWV data set. Distinctive patterns of enhanced autocorrelation are observable within the tropics and subtropics, in particular in the southern Pacific Ocean, with values reaching up to about 0.5. Towards higher latitudes, the distribution of the autocorrelation becomes spottier, and the values decrease to about 0.

After the calculation of the autocorrelation for each grid cell, the AR(1) model can be prepared via the transformation matrix P, as follows:

$$\begin{array}{}\text{(7)}& \mathbf{P}=\left[\begin{array}{ccccc}\sqrt{\mathrm{1}-{\mathit{\varphi }}^{\mathrm{2}}}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}& \mathrm{0}\\ -\mathit{\varphi }& \mathrm{1}& \mathrm{0}& \mathrm{⋮}& \mathrm{0}\\ \mathrm{0}& -\mathit{\varphi }& \mathrm{1}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{\ddots }& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{\cdots }& \mathrm{0}& -\mathit{\varphi }& \mathrm{1}\end{array}\right].\end{array}$$

For the case of the first element in the matrix, the AR(1) model cannot be constructed. Thus, the influence of the autocorrelation is approximated by $\sqrt{\mathrm{1}-{\mathit{\varphi }}^{\mathrm{2}}}$. If the time series has a gap between index t and t−1 (i.e. $X_t-X_{t-\mathrm{1}}>\mathrm{1}$), then the autocorrelation ϕ in Eq. (7) is set to 0 for this element.

Finally, the matrix P is then used to transform the fit function of Eq. (2) into the autocorrelation space as follows:

$$\begin{array}{}\text{(8)}& \mathbf{P}{Y}_t=Y_t^{\prime }=\mathbf{P}({\mathbf{M}}_{t}x+N_t)={\mathbf{M}}_{\mathbf{t}}^{\prime }x+{\mathit{\epsilon }}_t.\end{array}$$

The system of linear equations in Eq. (8) can then be solved by simple linear algebra in which the fit errors of the estimators already include the contribution from the autocorrelation of the noise.

One limitation of the AR model is the assumption of stationarity of the variance. Although this limitation can be overcome by using ARMA (autoregressive moving average) or ARIMA (autoregressive moving integrated moving average) processes, the determination and application of these models (for example, in the transformation of the linear equation system of the fit function) is highly nontrivial, especially for the case of unevenly spaced time series. Although an ARMA(1,1) process would be possible in the case that the lag 1 and lag 2 coefficients of the autocorrelation function have the same sign (e.g. Foster and Rahmstorf, 2011), this condition is not always given in our case. Thus, we have decided to stay with the AR(1) process.

3  Trend results

At this point, we would like to note that, especially in the high latitudes, complete temporal coverage within the MPIC OMI TCWV data set is not always given. For example, the winter months are often missing because no satellite measurements are available due to the seasonal solar cycle or ice cover. Thus, the trends shown are not representative for the entire year, but only for part of it, and should be interpreted with caution. However, we would still like to present the results, as these regions are of great interest in climate research. A map depicting the fractional temporal coverage is provided in Fig. S1 in the Supplement.

Moreover, when investigating climatological trends of TCWV on a local scale, these are also influenced by changes in atmospheric dynamics and should therefore be judged with caution. Nevertheless, they can still provide us with information about changes in the large-scale TCWV distribution.

3.1  OMI TCWV trends

To obtain reliable results, the trend analysis is performed only for grid cells whose time series cover at least half of the complete time period of interest. The results of the trend analysis of the OMI TCWV data set for the time range from January 2005 until December 2020 are illustrated in Fig. 2.

Figure 2 Global distributions of TCWV trends (2005–2020) derived from the MPIC OMI TCWV data set. Panels (a) and (b) depict the calculated absolute and relative TCWV trends, respectively. Panels (c) and (d) depict significant absolute and relative trends, respectively, after the application of the Z test. Panels (e) and (f) depict significant absolute and relative trends, respectively, after the application of the Z test and the false discovery rate (FDR) test. Grid cells for which no trend could be calculated (a, b) and/or for which the trends do not fulfil the significance criteria (c–f) are coloured grey.

The top row shows the absolute trends b (Fig. 2a) and the relative trends $\frac{b}{m}$ (Fig. 2b), respectively. Overall, increasing TCWV amounts are obtained. The absolute trends show high values in the equatorial Pacific and Southeast Asia, and the relative trends reveal high values in North America, the northern Pacific, and Southeast Asia. However, negative values in the TCWV trends can also be observed, e.g. in the region of the South Pacific convergence zone, southern Africa, Brazil, and the equatorial Atlantic. Altogether, we obtain a global area-weighted (i.e. weighted by the cosine of the latitude) mean absolute TCWV trend of +0.054 kg m⁻² yr⁻¹ and a relative TCWV trend of approximately +0.21 % yr⁻¹. We also obtained distinctively high trend values over mountains such as the Himalayas and Andes. However, these high values are likely artefacts due to uncertainties of the satellite retrieval, for example, in the input data for the ground elevation. Thus, we decided to filter these artefacts and only show grid cells for which the mean ground elevation is lower than 3000 m above mean sea level, based on the GMTED2010 elevation data set (Danielson and Gesch, 2011).

The linear least squares fit assumes that errors of the estimators are normally distributed. Thus, we can perform a Z test from the fit results and determine which trends are statistically significant or not. For our purposes, we choose a significance level of 5 %, for which the Z test requires that |b|≥1.96σ_b (see Fig. 2c and d). Furthermore, to account for test multiplicity and field significance, we additionally perform a false discovery rate (FDR) test (Benjamini and Hochberg, 1995; Wilks, 2006, 2016). Because the OMI TCWV data set also shows a high spatial autocorrelation (see Appendix B), we follow the recommendations in Wilks (2016) and choose a significance level of 2.5 % for the FDR test.

The remaining trends are given in the bottom row of Fig. 2, with the absolute and relative trends in panels (e) and (f), respectively. From about 12 500 trends originally classified as significant according to the Z test, approximately 4700 grid cells still remain significant after the application of the FDR test, and almost all of them reveal a positive TCWV trend, in particular over the Pacific Ocean, East Asia, and parts of the U.S. East Coast.

In addition to the TCWV trends, we also analyse the trends of the individual components of the DOAS retrieval, i.e. the slant column density (SCD) and the air mass factor (AMF), where TCWV = SCD/AMF. These additional analyses reveal that the TCWV trends are mainly determined by trends in the SCD, i.e. by increasing or decreasing H₂O absorption due to changing atmospheric water vapour content, respectively. The trends of the inverse AMF (i.e. 1/AMF) are generally negative but also distinctively weaker (about 3–4 times) than the SCD trends and thus have only a moderate influence on the overall TCWV trends. More details on these analyses are given in Appendix C.

Figure 3 Global distributions of TCWV trends (2005–2020) derived from the MPIC OMI TCWV data set. Panels (a) and (b) depict the calculated relative TCWV trends with and without teleconnection indices, respectively, in the trend analyses. Panel (c) depicts the differences between the trend results with teleconnections in the analysis minus the trend results without teleconnections in the analysis. Grid cells for which no trend could be calculated are coloured grey.

To highlight the influence of teleconnections on the trend results for the OMI TCWV data set, we also perform the trend analysis not accounting for them. The resulting trends and their difference are shown in Fig. 3. While overall the spatial distributions of the relative trends (Fig. 3a and b) look quite similar, distinct patterns emerge when looking at the trend difference (Fig. 3c). For instance, the typical PMM and ENSO teleconnection patterns are clearly visible (e.g. dipole structure over the maritime continent in the case of ENSO). Consequently, the resulting deviations are particularly strong in the tropical and subtropical Pacific and can reach values as high as the relative trends themselves. We have also tested other AR models with lag =2, 3, 6, and 12 and found that the trend results and the distributions of the significant trends differ only slightly from those using an AR(1) model. The corresponding trend results can be found in Fig. S2 in the Supplement.

3.2  Intercomparison to trends of other TCWV data sets

To verify the OMI TCWV trends and to detect potential shortcomings within the OMI TCWV data set, we performed the analyses also for monthly mean TCWV data from the reanalysis model ERA5 (Hersbach et al., 2019, 2020). For this purpose, the ERA5 TCWV data set is gridded on a ${\mathrm{1}}^{\circ }\times {\mathrm{1}}^{\circ }$ lattice. Moreover, to account for OMI's observation time (13:30 LT), we only take into account ERA5 monthly mean values between 13:00–14:00 LT.

Figure 4 Global distributions of relative TCWV trends derived from the OMI TCWV data set (a, c) and ERA5 (b, d). Panels (a) and (b) depict all calculated relative TCWV trends, and panels (c) and (d) show the corresponding significant trends remaining after the application of the Z test and FDR test. Grid cells for which no valid trend could be calculated are coloured grey.

The resulting maps of the relative trends are given in Fig. 4. Overall, the trend results of OMI and ERA5 agree well to each other, as both all and only significant relative trend results (top and bottom rows in Fig. 4, respectively) have similar strengths and also show similar global distributions. Nevertheless, the OMI TCWV trends reveal slightly stronger increases over parts of East Asia (which are also classified as significant) and South America and are in general less smooth than the ERA5 results. Similar findings can be obtained for the absolute trends, which are available in Fig. S3 in the Supplement.

Figure 5 Global distributions of relative TCWV trends derived from the OMI TCWV data set (a, c) and GOME-Evolution (b, d) for the time range from January 2005 to December 2015. Panels (a) and (b) depict all calculated relative TCWV trends, and panels (c) and (d) show the corresponding significant trends remaining after the application of the Z test and FDR test. Grid cells for which no valid trend could be calculated are coloured grey.

In addition to ERA5, we also compare the trend results to trends from the TCWV satellite product GOME-Evolution (Beirle et al., 2018). Since the GOME-Evolution product is only available until 2015, we modified the time range accordingly, i.e. the results for the relative trends shown in Fig. 5 (and for the absolute trends in Fig. S4) correspond to a time range from January 2005 to December 2015. While the distributions of the relative trends have quite similar patterns and partly similar magnitudes, striking differences can be seen in some regions. For example, the OMI trends in the tropical Pacific North America or the Arabian Peninsula are much higher than the GOME-Evolution trends. Also, overall, many more trends are classified as significant for OMI than for GOME-Evolution.

Nevertheless, considering that the GOME-Evolution product retrieves total column water vapour in the visible red spectral range, uses a different vertical column density (VCD) conversion scheme (see also Wagner et al., 2003, 2007; Grossi et al., 2015), and observes the atmosphere at an earlier overpass time (around 10:00 LT), the good agreement in the trend results further confirms the reliability of the findings of the OMI TCWV trend analysis.

Furthermore, we made additional comparisons to the results of past studies. From these comparisons, several differences in the strength and spatial distribution of TCWV trends emerge. The reasons for these differences are, on the one hand, the consideration of different time periods and, on the other hand, also different methods of analysis. Further details about these comparisons can be found in Appendix D.

4  Trends in relative humidity

In this section, we investigate to what extent the assumption of constant relative humidity is given at the local scale. For this purpose, we make the following assumptions. First, we assume that the relative changes in TCWV correspond to those in near-surface specific humidity q_s, i.e. $\frac{\mathrm{dTCWV}}{\mathrm{TCWV}}\approx \frac{\mathrm{d}{q}_{\mathrm{s}}}{q_{\mathrm{s}}}$. This assumption should be fulfilled, since TCWV is directly connected to the specific humidity via its vertical integral, and approximately 60 % of the TCWV is located within the planetary boundary layer. Second, we also assume that relative changes in specific humidity correspond to changes in water vapour pressure, i.e. $\frac{\mathrm{d}q}{q}\approx \frac{\mathrm{d}e}{e}$ (assuming that relative changes in surface air pressure are negligible, i.e. $\frac{\mathrm{d}{p}_{\mathrm{s}}}{p_{\mathrm{s}}}\ll \frac{\mathrm{d}e}{e}$). Given the aforementioned assumptions and that the water vapour pressure e can be described as $e=\mathrm{RH}\cdot E$, we can derive the relative changes in relative humidity (RH) by combining the relative TCWV trends with trends in surface air temperature T, as follows:

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\mathrm{d}{q}_s}{q_{\mathrm{s}}}}\approx {\displaystyle \frac{\mathrm{d}e}{e}}& {\displaystyle }={\displaystyle \frac{\mathrm{dRH}}{\mathrm{RH}}}+{\displaystyle \frac{\mathrm{d}E}{E}}\text{(10)}& {\displaystyle }\to {\displaystyle \frac{\mathrm{dRH}}{\mathrm{RH}}}& {\displaystyle }={\displaystyle \frac{\mathrm{d}{q}_s}{q_{\mathrm{s}}}}-{\displaystyle \frac{L_{\mathrm{v}}\left(T\right)}{R_{\mathrm{v}}}}{\displaystyle \frac{\mathrm{d}T}{T^{\mathrm{2}}}}\approx {\displaystyle \frac{\mathrm{dTCWV}}{\mathrm{TCWV}}}-{\displaystyle \frac{L_{\mathrm{v}}\left(T\right)}{R_{\mathrm{v}}}}{\displaystyle \frac{\mathrm{d}T}{T^{\mathrm{2}}}}.\end{array}$$

Thus, if RH is 50 %, a relative increase of 1 % indicates an absolute RH increase of 0.5 %. However, it should be noted that the largest uncertainties lie in the first assumption, i.e. slight under- or overestimations of the actual relative q_s changes will cause corresponding deviations in the relative RH changes.

Figure 6 Relative trends in relative humidity (RH) derived from the relative TCWV trends and the temperature trends from OMI and Berkeley Earth (a), from ERA5 (b), and from the data set HadISDH (Hadley Centre Integrated Surface Dataset of Humidity) (c) for the time range from January 2005 to December 2020. Grid cells for which no trend has been calculated are coloured grey.

Figure 6 depicts the resulting relative RH trends derived from the OMI TCWV trends in combination with the temperature trends from the Berkeley Earth temperature data record (Rohde and Hausfather, 2020), from ERA5, and from the relative RH trends from the HadISDH (Hadley Centre Integrated Surface Dataset of Humidity) surface relative humidity data set (Willett et al., 2014, 2020). In general, the results for OMI and ERA5 reveal a global (relative) increase in RH, in particular that the trends over ocean are widely positive. However, in all three data sets, distinctive decreasing trends are observable over land, for instance over Russia or southern Africa. Considering the differences in the selected time period and measurement source, the RH trends from OMI over land surface coincide well with the results from Dunn et al. (2017). The reduction in relative RH over land is likely related to a marked land–ocean contrast in warming (Simmons et al., 2010; Fasullo, 2012), besides various local factors such as changes in vegetation cover (Simmons et al., 2010). Over ocean, due to the direct link with sea surface temperature, the water vapour content can increase adequately to keep RH constant. Over land, this is usually only possible with a delay due to limited water availability, as water must first be transported there from ocean. Since the temperature also increases much more over land than over ocean, the decrease in RH might be due to the lack of an increased water supply from the ocean (Simmons et al., 2010).

Interestingly, we also find distinctive increases in RH in arid regions (e.g. over the Sahara) and distinctive decreases in humid regions (e.g. the tropical Pacific Ocean) within the OMI and the ERA5 results. Recently, Bourdin et al. (2021) investigated RH trends from the reanalysis models ERA5 and JRA-55 (Japanese 55-year Reanalysis) over the past 40 years and also found significant negative trends in the tropical lower troposphere.

Several studies have shown that global warming will lead to a further drying of dry regions (e.g. Sherwood and Fu, 2014), and wet regions will become even wetter (e.g. Held and Soden, 2006; Chou et al., 2013; Allan et al., 2010), leading to the simple paradigm of “dry gets drier, wet gets wetter” (DDWW; Chou et al., 2009). In addition, other studies show that changes in precipitation correlate very well with changes in ocean salinity, suggesting a “fresh gets fresher, salty gets saltier” pattern (Cheng et al., 2020, and references therein). Though most of these studies focus on changes in precipitation, our results for RH support the findings from Greve et al. (2014) and Byrne and O'Gorman (2018) in that the DDWW paradigm is not always fulfilled over land. Surprisingly, according to our results, this paradigm is not fulfilled even over the tropical Pacific Ocean, the region on which most of the concepts of the studies are based (e.g. Held and Soden, 2006). However, we would like to stress here that the time period studied is probably too short to question the paradigm.

5  Relationship between TCWV and precipitation

According to Bretherton et al. (2004) and Rushley et al. (2018), a nonlinear relationship between TCWV (or column relative humidity, respectively) and precipitation exists for the tropical ocean. Thus, given the TCWV and RH trend results, we expect to observe a decline or negative trend, in particular over the Pacific Ocean, along the Intertropical Convergence Zone. For the analysis of trends in precipitation, we use the monthly mean rain rates from the GPCP Version 3.2 Satellite-Gauge (SG) Combined Precipitation Data Set (Huffman et al., 2020). For the sake of consistency, we grid the GPCP data from a resolution of ${\mathrm{0.5}}^{\circ }\times {\mathrm{0.5}}^{\circ }$ to a ${\mathrm{1}}^{\circ }\times {\mathrm{1}}^{\circ }$ lattice.

Although precipitation climate data records (CDRs) allow a global analysis, they are subject to large uncertainties, as satellite and rain gauge observations do not have good spatiotemporal coverage, weak and short rain events are not well detected or even missed, and satellite retrievals can determine the rain rate only indirectly. Thus, deviations of about 50 % in the daily rain rate can occur, compared to in situ measurements (e.g. Prat et al., 2021). Nevertheless, Prat et al. (2021) show that, over accumulation periods of month or years, precipitation CDRs perform satisfactorily. Moreover, Prat et al. (2021) used an older GPCP version (v2) than ours in their evaluation study.

Figure 7 Global distribution of relative RH trends derived from the OMI TCWV data set (time range from 2005 to 2020; a; same as in Fig. 6a) and of trends in precipitation derived from the GPCP v3.2 monthly mean data set for the time range from January 2005 to December 2020. Panel (b) depicts all rain rate trends, and panel (c) shows only those that are considered significant after applying the Z test (to the significance level of 5 %) and a FDR test (see also Appendix B). Grid cells for which no valid trends have been calculated are coloured grey.

Figure 7 depicts the obtained trends in precipitation and the relative RH trends from OMI. Comparing the trend distributions of the monthly mean rain rates to the relative RH trends, negative and positive trends in precipitation and RH match quite well over the tropical and subtropical ocean, especially over the tropical Pacific and the northern subtropical Atlantic. While over land within the subtropics an acceptable match can be determined in some regions (e.g. southern Africa and Brazil), the patterns of the relative RH and rain rate trends no longer match well towards mid and high latitudes (e.g. in North America), likely because, in these regions, the rain rate is mainly determined by atmospheric dynamics (cyclone or storm tracks) rather than thermodynamics. Furthermore, the distinctive relative increases in RH in the mountainous regions of South America (Andes) and northern India (Himalayas) are likely due to the inadequacies in the OMI TCWV satellite data caused by the complex topography (see also Sect. 3.1).

Trenberth (2011) and Trenberth and Shea (2005) analysed local correlations between precipitation and surface temperature for cold and warm seasons and reported mainly positive correlations over ocean and negative correlations year round over land throughout the tropics. However, over ocean, the correlations also depend on whether the (sea) surface temperature is driven by the ocean or by the atmosphere (Trenberth and Shea, 2005). While in some regions of the subtropics we can also find this high correlation in the trend patterns of precipitation and surface temperature (e.g. increase in the precipitation in the northern subtropics in the eastern Pacific or decrease in the subtropical Atlantic over ocean; decrease in Brazil or southern Africa over land), we cannot find a direct link for the striking negative precipitation trends in the equatorial Pacific. However, it should also be taken into account that a large part of the precipitation trends are not statistically significant.

Overall, the discrepancies between our observations and the expected changes in the hydrological cycle show that accurate observations and long-term monitoring of the Earth's hydrological cycle and atmosphere on the global scale from multiple remote sensing and in situ platforms are essential to clarify this important aspect.

6  Changes in the atmospheric water vapour residence time

Another key diagnostic of the hydrological cycle is the atmospheric water vapour residence time (WVRT). The WVRT can contribute to a better understanding of changes in dynamic and thermodynamic processes within a changing climate (Trenberth, 1998; Gimeno et al., 2021). For instance, an increase in WVRT suggests that the length of the atmospheric moisture transport increases, i.e. the distance between moisture sink and source regions (Singh et al., 2016). Several different metrics exist for quantifying the WVRT (van der Ent and Tuinenburg, 2017; Gimeno et al., 2021), bearing in mind that the WVRT distribution or the lifetime distribution (LTD) is exponential on the local scale, and thus, the mean value is strongly influenced by a few high values (van der Ent and Tuinenburg, 2017; Sodemann, 2020). Ideally, one would determine the LTD for each grid cell for each month from backward trajectories and then examine their changes or trends. However, this would be well beyond the scope of this paper.

Thus, for our purpose, and for the sake of simplicity, we focus on the so-called depletion time constant (DTC) and the turnover time (TUT). The TUT describes the global average mean age of precipitation and can be calculated as the ratio of TCWV to precipitation P, as follows:

$$\begin{array}{}\text{(11)}& \text{TUT}={\displaystyle \frac{\stackrel{\mathrm{‾}}{\text{TCWV}}}{\stackrel{\mathrm{‾}}{P}}},\end{array}$$

where the bar indicates global average. Typically, the TUT varies between values of 8 to 10 d and is expected to increase by 3 % K⁻¹–6 % K⁻¹ (Gimeno et al., 2021, and references therein). Analogously, the DTC is defined as the local ratio of TCWV to precipitation, as follows (e.g. Trenberth, 1998):

$$\begin{array}{}\text{(12)}& \text{DTC}={\displaystyle \frac{\text{TCWV}}{P}}.\end{array}$$

The DTC values might vary substantially from TUT, but the global precipitation weighted average is equal to TUT (Gimeno et al., 2021).

For our investigations of trends in DTC, we combine the regridded GPCP data set from Sect. 4 and the OMI TCWV data set and perform the trend analysis scheme from Sect. 2.2 to the monthly DTC values for the time range 2005 to 2020. To ensure numerical stability, we only consider monthly rain rates greater than 0.25 mm d⁻¹. As a result, large parts of the subtropical oceans and deserts are excluded from the analysis.

The results of the DTC trend analyses are depicted in Fig. 8. On average, we typically obtain mean DTC values between 5–10 d in the areas where rain occurs (Fig. 8a). In the subtropical dry zones, values of around 30 d and well above are found. In terms of absolute DTC trends, the most striking patterns are in the northern subtropical Atlantic, with strong increases, and in the northern subtropical western Pacific, with strong decreases. In comparison, the distribution of relative DTC trends is much spottier, but overall, in addition to the patterns already mentioned, we obtain distinctive increases in DTC on the U.S. West Coast, in Europe, Russia, and in the eastern Pacific.

For our investigations of trends in TUT, we first calculate global averages of the regridded GPCP data set from Sect. 4 and the OMI and ERA5 TCWV data sets between 60^∘ S and 60^∘ N for each month, then combine the time series of global averages, and, finally, perform the trend analysis for the TUT time series for the time range from 2005 to 2020.

Altogether, we find an increase in the global TUT for OMI and ERA5 of approximately +0.02 d yr⁻¹, with TUT mean values of around 9.7 d for OMI and 8.8 d for ERA5. Combining the long-term relative trends in TUT and trends in surface air temperature, we can estimate the sensitivity of TUT to global warming r, i.e. $r=\frac{\mathrm{\Delta }\mathrm{TUT}}{\mathrm{TUT}}/\mathrm{\Delta }T$. For the case of OMI and Berkeley Earth, we find a TUT sensitivity of around 8.4 % K⁻¹ and for ERA5 of around 8.8 % K⁻¹, which is higher than the results of 3 % K⁻¹–6 % K⁻¹ pooled in Gimeno et al. (2021).

Figure 8 Global distribution of DTC trends for the time range of January 2005 to December 2020. Panel (a) depicts the distribution of the mean DTC. Panels (b) and (c) depict the absolute and relative DTC trends. Grid cells for which no valid trend has been calculated are coloured grey.

7  Summary

In this study, we analysed global trends within a long-term data set of total column water vapour (TCWV) retrieved from multiple years of OMI observations for the time period January 2005 until December 2020 and considered the effects of autocorrelation of the residuals within the analysis scheme. The results of the analyses were then put into context with trends from additional TCWV data sets, like from the GOME-Evolution project or from the reanalysis model ERA5, and overall very good agreement was found. In a next step, based on the relative OMI TCWV trends, trends in relative humidity were derived and put into the context of the assumption of invariant relative humidity. Moreover, under consideration of the relationship between (column) relative humidity and precipitation, the patterns of the relative RH trends have been compared to rain rate trends. Also, the changes in the water vapour residence time and its response to changes in surface air temperature were investigated.

The trend analysis reveals an increase in TCWV of approximately +0.054 kg m⁻² K⁻¹ or +0.21 % yr⁻¹ globally for the time period of January 2005 until the end of 2020. To determine if trends are significant or not, a Z test and a false discovery rate test are applied to the trend results. After application of these significance criteria, almost all remaining trends are positive and distributed across the globe. However, particular spatial patterns remain, for instance within the region of the northern subtropics of the eastern Pacific. Overall, the relative OMI TCWV trends agree well to the corresponding trends from ERA5 and from the GOME-Evolution data set.

To analyse if the assumption of temporally invariant relative humidity is fulfilled on the local scale, we derived relative trends in relative humidity (RH) from the TCWV trends. All in all, we obtain that RH increases distinctively over large areas of the ocean and land surface. However, over both surface types relative decreases can also be well identified in some areas. Interestingly, relative decreases and increases in RH are not limited to arid and humid regions, respectively. For instance, our analysis reveals relative increases in RH over the (arid) Sahara desert and decreases in RH over the (humid) tropical Pacific Ocean. Within the tropics, we also find that the patterns of decreasing RH trends match those of decreasing precipitation quite well, especially within the tropical Pacific Ocean.

Combining the TCWV and rain rate data sets, changes in the water vapour residence time (WVRT) have been investigated. Overall, an increase in the turnover time of about 0.02 d yr⁻¹ has been observed. Together with the long-term trends in surface temperature, we estimate a TUT sensitivity to global warming of around 8.4 % K⁻¹, which is 2 to 3 times higher than the values provided in Gimeno et al. (2021).

All in all, our results show that several challenges still remain for a better understanding of the atmospheric hydrological cycle and even new questions arise regarding the complex interactions between air temperature, water vapour, precipitation, and atmospheric dynamics. The differences between observed and expected changes in the hydrological cycle show that simplified assumptions are not always valid (e.g. invariant relative humidity). Also, our observed, much higher global sensitivities of individual parameters of the hydrological cycle (i.e. TUT) to changes in surface temperature raise the question of what effects can be expected at the local scale (e.g. precipitation) with further increasing temperatures, especially with regard to changes in the global circulation such as the expansion of the Hadley cell towards higher latitudes (e.g. Staten et al., 2018).

With regard to TCWV retrievals in the visible blue spectral range, there is great potential for extending the OMI TCWV data set with further satellite data (e.g. from TROPOMI or GOME-2) and combining it with future missions from geostationary satellites, such as GEMS or Sentinel-4, which will also allow for investigations of (semi-)diurnal TCWV cycles.

Appendix A:  Influence of the autocorrelation on the trend results

To address the influence of the autocorrelation on the trend results for the OMI TCWV data set, we perform the trend analysis not accounting for it. The panels in Fig. A1 illustrate the difference in the absolute (Fig. A1a) and relative (Fig. A1b) trend results (i.e. the difference of the results with accounting minus the results without accounting for the influence of the temporal autocorrelation). For high and mid latitudes, the differences are close to zero, indicating that the influence of the autocorrelation on the trend results is negligible. However, within the subtropics and tropics, distinctive deviations are observable, especially in the regions where the autocorrelation is high (e.g. the Pacific Ocean; see also Fig. 1). For the case of the relative trends (Fig. A1b), the deviations can reach up to 0.05 % yr⁻¹ (which is around 10 % of the magnitude of the relative trends in the affected regions) and consequently can cause wrong signs in the trend estimation (i.e. indicating a negative instead of a positive trend).

Figure A1 Difference between trends of the MPIC OMI TCWV data set (2005–2020) with accounting minus not accounting for the influence of autocorrelation, where panel  (a) shows the absolute trends and panel (b) the relative trends.

Appendix B:  Spatial autocorrelation within the MPIC OMI TCWV and GPCP data set

The significance level at which the false discovery rate test method in Sect. 3.1 is performed depends on the degree of spatial autocorrelation. Thus, for every timestamp within the MPIC OMI TCWV data set, the spatial autocorrelation is calculated from the global TCWV distribution for grid point separations up to 7000 km.

Figure B1 Spatial autocorrelation as function of the great circle distance of the MPIC OMI TCWV (a) and the GPCP data set (b). The black dots represent the results of the analysis of the spatial distribution for each time step in the respective data set. The solid red lines illustrate the fit result of $f\left(x\right)=e^{-a|x{|}^b}$.

Figure B1a illustrates the spatial autocorrelation of the OMI TCWV data set as a function of grid point separation. The red solid line is the fit result of $f\left(x\right)=e^{-a|x{|}^b}$, with the grid point separation distance x. For the OMI TCWV data set, we calculated a value of a≈0.098 and b≈1.88, which equals an e-folding distance of approximately 3.43×10³ km. According to Wilks (2016), this e-folding distance indicates a strong spatial dependency. Consequently, we follow the recommendations of Wilks (2016) and set, for the FDR test, the significance level to 2.5 % instead of 5.0 %.

The same procedure was applied to the GPCP data set, the results of which are shown in Fig. B1b. For the fit results, we obtain a≈0.58 and b≈1.28, which correspond to an e-folding distance of 1.53×10³ km, which is thus less than half as large as that of TCWV. Accordingly, the spatial dependence is not so strong, and the significance level for the FDR test can remain at 5 % for the GPCP data set.

Appendix C:  Trends of individual retrieval parameters

Here, we investigate to what extent the relative TCWV trends are due to geophysical changes in the water vapour content or due to changes in the retrieval input parameters. For DOAS retrievals, the TCWV amount is derived via the quotient of the integrated concentration along the light path (so-called slant column density, SCD) and the so-called air mass factor, AMF, i.e. TCWV = SCD/AMF. Thus, the relative trends of these two quantities were calculated following the analysis scheme in Sect. 2.2. For the case of the SCD, we use the geometrical VCD (vertical column density), which is simply the SCD divided by the geometrical air mass factor (which remains constant over time).

Figure C1 Global distributions of relative trends of the TCWV (a), geometrical vertical column density (VCDgeo; b), and the inverse of the air mass factor (1/AMF; c) for the time period from January 2005 to December 2020. Grid cells for which no trend has been calculated are coloured grey.

The global distributions of the relative trends of both quantities are illustrated in Fig. C1b and c and the relative TCWV trends (in Fig. C1a). The distribution and strength of the geometrical VCD (Fig. C1b) largely coincide with the distribution of the relative TCWV trends (Fig. C1a). The trends of the inverse AMF (1/AMF, Fig. C1c), on the other hand, are, in general, much weaker than the SCD trends (approx. 3–4 times weaker) and do not follow the TCWV trend distribution. However, it occasionally happens that the relative inverse AMF trends either weaken or cancel the SCD trends (e.g. North America or northeastern Asia) or even strengthen them (e.g. around the Arabian Peninsula). Overall, we conclude that the relative TCWV trends are mainly determined by the SCD trends, which consequently means that TCWV trends are mainly due to an increase in atmospheric water vapour concentration.

Figure C2 Absolute trends of the retrieval input parameters for the calculation of the air mass factor for the time period from January 2005 to December 2020. (a) Effective cloud fraction. (b) Cloud-top height. (c) Surface albedo. Grid cells for which no trend has been calculated are coloured grey.

In addition to the trends of the SCD and AMF, we also analyse the trends of the AMF input parameters, i.e. the effective cloud fraction (CF), the cloud-top height (CTH), and the surface albedo. The corresponding global distributions are depicted in Fig. C2. Here, it is important to mention that the MPIC OMI TCWV data set only includes mostly clear-sky observations (i.e. CF <20 %), so the calculated trends of the cloud input parameters are very likely not representative for the actual cloud trends of the atmosphere. For CF (Fig. C2a), we obtain, in general, decreasing trends around −0.1 % yr⁻¹ globally, except for the Indian subcontinent and some individual locations. For the input CTH (Fig. C2b), no clear trend pattern is observable, except for slight increasing trends over the tropical landmasses with values around +0.03 km yr⁻¹. As expected for the surface albedo (Fig. C2c), no trends are observable over ocean, as a static monthly albedo map has been used here. Over land, however, strong varying trends can be found in the high latitudes of the Northern Hemisphere, with absolute values higher than 0.2 % yr⁻¹. Nevertheless, these strong albedo trends in the Northern Hemisphere are typically not significant.

Appendix D:  Intercomparison to trends from other studies

In the following, we compare our results of relative TCWV trends for the time range 2005–2020 to trends presented in previous studies and investigate which TCWV trends are significant within the respective time range of the previous studies. It is particularly important to note that TCWV trends from different time periods have been investigated. For the sake of completeness, the global distributions of the absolute trends for the same data sets and time ranges are available in Fig. S5 in the Supplement.

Figure D1 Global distributions of relative TCWV trends of OMI (2005–2020; a, b) and ERA5 for the following different time periods: (c, d) 1988–2003, (e, f) 1996–2006, and (g, h) 1995–2011. Panels in the left column illustrate all calculated trends, and panels in the right column illustrate statistically significant trends after the application of a Z test and a FDR test. Grid cells for which no valid trend has been calculated are coloured grey.

Trenberth et al. (2005) analysed trends from the RSS SSM/I data for the time period of 1988 to 2003. While the patterns generally match quite well, the trends often have opposite signs. In our period (2005–2020), the trends are mainly positive, whereas in the period of this study (1988–2003), the trends are mainly negative. This is particularly visible in the eastern Pacific. However, we were unable to identify any significant trends for this period (see Fig. D1d). Overall, however, the trends of Trenberth et al. (2005) are in very good agreement with the trends we have determined for this period (compare Fig. 11 in their paper).

Mieruch et al. (2008) investigated TCWV trends from 1996 to 2006, using a TCWV data set created from measurements of GOME and SCIAMACHY, using the AMC-DOAS method (Air Mass Corrected Differential Optical Absorption Spectroscopy; Noël et al., 2004). In contrast to the comparison with Trenberth et al. (2005), almost no similarities are discernible either in the spatial patterns or in the strength of the trends. Overall, the spatial distribution is not as smooth as in the other periods studied and is distinctively spottier. This is probably due to the fact that the period studied is quite short and that there was also a strong El Niño event in 1997/1998. Compared to the results in Mieruch et al. (2008, Fig. 5 in their paper), our results for the same period find only few similarities, even for the case of the significant trends. For example, the trends of Mieruch et al. (2008) are sometimes 4 to 6 times higher than ours for the same period.

More recently, Wang et al. (2016) also investigated TCWV trends for the time period from 1995 to 2011 for a TCWV data set combining measurements from radiosondes, GPS radio occultation, and microwave satellite instruments. As for the comparison to Trenberth et al. (2005), our findings and the findings from Wang et al. (2016) share many similarities but also several discrepancies. Wang et al. (2016) find a “sandwich” shape in the tropical and subtropical Pacific, with positive trends in the region of the Intertropical Convergence Zone bounded by two bands of negative trends. In contrast, the OMI TCWV trends also suggest a sandwich shape but with opposite signs to Wang et al. (2016), i.e. negative trends bounded by positive trends. Such opposite findings also occur over parts of the Indian subcontinent, the Arabian Peninsula, and South America. However, for central Europe and parts of Asia, good agreement for the trend patterns is found.

For the comparisons of our results to the findings of Trenberth et al. (2005), Mieruch et al. (2008), and Wang et al. (2016), one explanation for the differences may be the different time periods of investigations (1988 to 2003, 1996 to 2006, and 1995 to 2011 vs. 2005 to 2020). Figure D1c–h illustrate the relative TCWV trends derived from the ERA5 data set for the aforementioned time periods. Although only the time periods have been changed, clear differences can indeed be identified in both the distribution and the strength of the trends. Furthermore, these trend distributions agree very well with the results of the three previously mentioned studies. Nevertheless, different methodologies of observations or different methods for the trend calculation may also be a cause for the discrepancies. For instance, we explicitly account for the influence of ENSO by including the ONI and TNI index into our analysis scheme (see also Sect. 2.2 and Sect. 3.1), whereas Mieruch et al. (2008) explicitly filtered the time around the strongest ENSO signal.

Combining that the detected trends for ERA5 and the GOME-Evolution data set agree well to the findings from the OMI TCWV data set (see Sect. 3.2) but that the comparisons to the results from other trend analysis studies show systematic differences, it is evident to not only compare trends for the same time periods but also to ensure that the same methodology for the trend analysis is used. As a lot of different methods exist for estimating trends in environmental data sets, it would be particularly interesting to evaluate which trend analysis scheme performs best and should be recommended for future studies. However, such an evaluation study is beyond the scope of this paper.