Tropical rainfall subseasonal-to-seasonal predictability types

Tropical rainfall is mostly convective and its subseasonal-to-seasonal (S2S) prediction remains challenging. We show that state-of-art model forecast skill 3 + 4 weeks ahead is systematically lower over land than ocean, which is matched by a similar land-ocean contrast in the spatial scales of observed biweekly rainfall anomalies. Regional differences in predictability are then interpreted using observed characteristics of daily rainfall (wet-patch size, mean intensity as well as the strength of local S2S modes of rainfall variation), and classified into six S2S predictability types. Both forecast skill and spatial scales are reduced over the continents, either because daily rainfall patches are small and poorly organized by S2S modes of variation (as over equatorial and northern tropical Africa), or where the daily mean intensity is very high (as over South and SE Asia). Forecast skill and spatial scales are largest where daily rainfall is synchronized by intraseasonal (such as the Madden-Julian Oscillation) as well as interannual ocean-atmosphere modes of variation (such as El Niño-Southern Oscillation), especially over northern Australia and parts of the Maritime Continent, and over parts of eastern, southern Africa and northeast South America. The oceans exhibit the highest skill and largest spatial scales, especially where interannual (central equatorial Pacific) or intraseasonal (central and eastern Tropical Indian Ocean and Western Pacific) variability is largest. These results provide a relevant regional typology of the potential drivers and controls on S2S predictability of tropical rainfall, informing intrinsic limits and possible improvements toward useful S2S climate prediction at regional scale.


INTRODUCTION
Tropical rainfall variability strongly impacts local societies, especially through the annual rhythm of the monsoons with distinct wet and dry seasons, which often undergo pronounced interannual and intraseasonal fluctuations. 1,2 A major goal of current international research efforts is to improve forecasts and understanding on the subseasonal to seasonal (S2S) range, 2 weeks to a season ahead, filling the gap between medium range weather forecasts and seasonal forecasts. [3][4][5][6][7] This timescale includes the meridional seasonal migration of the intertropical convergence zone (ITCZ) and the monsoons, together with the primary modes of tropical S2S variability, especially the intraseasonal Madden Julian Oscillation (MJO [8][9][10], and El Niño-Southern Oscillation (ENSO 11,12 ) on the seasonal to interannual scales. [11][12][13][14] However, S2S prediction of tropical rainfall over land remains a major challenge, and the geographical distribution of its predictability and the factors that limit it are poorly understood. 6,7,15,16 In particular, the S2S scale encompasses interactions between intrinsic monsoon dynamics and the MJO leading to monsoon intraseasonal oscillations typified by active and break phases and changes in the onset date. [17][18][19] The famous Koppen-Geiger classification of local mean climates was developed over 100 years ago 20 and is still in widespread use in atlases and beyond. 21 However, an analogous typology for rainfall variability and predictability is still lacking. The goal of this paper is to make a relevant geographical classification of tropical rainfall predictability and its controlling factors from both S2S forecast models and observed data, thereby identifying the countries and regions that could benefit most from S2S prediction through development of effective weather and climate services. 22 Climate prediction is a signal-to-noise (SN) problem where the signal (noise) is the predictable (unpredictable) fraction of the climate variation at a given timescale. [23][24][25] S2S predictions of weekly or biweekly average rainfall 5,6 are possible where the signal lent by predictable modes of climate (especially MJO and ENSO) exceeds the unpredictable noise, largely associated with sub-weekly weather, [26][27][28][29] provided that the forecast model can capture the predictable signals. This competition between signal and noise can be estimated from S2S forecast model outputs alone, and by confronting the model output with observations over a reforecast period (model skill). On the other hand, S2S signals and noise can also be assessed empirically, purely from observed rainfall data. Larger spatial scales of biweekly average observed rainfall anomalies imply larger S2S predictability because the predictable climate drivers (MJO and ENSO) are large scale and persistent at this timescale, while the unpredictable weather noise in tropical convective rainfall is small scale and random in time and space. In this way, the spatial scales of rainfall time averages can empirically measure the upper bound of S2S predictability at local or gridbox scale from observed data alone. 30,31 Determining the potential for S2S prediction at local to regional scales requires mapping both signals and noise, as well as their underlying physical controls. We combine both model-based and observations-based estimates and focus here on biweekly rainfall amounts at a 2-4 week forecast lead time (15-28 days ahead, also referred to as forecast weeks 3 + 4). We first assess skill and ensemble spread from a 20-year ensemble of 11 runs of ECMWF model reforecasts. These model-based estimates of predictability are then interpreted geographically by clustering various empirical estimates of signal and noise, based on observed rainfall data. S2S predictable signals are defined empirically using three measures: the fraction of biweekly rainfall variance in the 7-20 day and 20-90 day intraseasonal, and interannual frequency bands. The noise is expressed empirically using four measures: the average area or "patch-size" of daily rainfall events, the mean daily rainfall intensity, spatial autocorrelation of biweekly rainfall anomalies, and its subweekly variance fraction.

RESULTS
Model skill and predictability and observed spatial scales both peak over the equatorial oceans ECMWF model forecast skill-defined here by the Spearman correlation between the ensemble mean of the 11 runs and the observations-and potential predictability 32 (PP)-defined here by the averaged Spearman correlation between each run vs the ensemble mean-for week 3 + 4 rainfall are both highest over the equatorial Pacific ENSO core region, and are smaller over land than over the ocean 5,6 ( Fig. 1). The pattern correlation (PC hereafter) between both maps equals 0.91. Over land, a relatively large PP is not always associated with high skill, such as over western India or Southern Africa, which may reflect poor simulation of predictable regional circulation patterns or a relatively large amplitude of internal atmospheric variations due to chaotic dynamics. 33,34 There is a striking similarity between the geographical pattern of model skill (Fig. 1a) and spatial scales (Fig. 2) since the largest spatial co-variations in skill and spatial scale (PC = 0.7) are between land (relatively low skill and spatial scale) and ocean (relatively high skill and spatial scale). Spatial scales are significantly smaller across the landmasses, especially over northern and equatorial Africa, southern Asia between the Gangetic plain and Indochina, western South America and most of central America. Scales are larger over the equatorial Pacific and around the Maritime Continent. Pockets of larger spatial scales over eastern South America or eastern Africa (Fig. 2) usually coincide with higher skill in Fig. 1a. These results demonstrate that observed spatial scales of biweekly rainfall variability are usually larger where model forecast skill (and PP) are large and vice versa, with oceanic rainfall variations being more predictable and larger scale, than over land.
What explains these skill/scale geographical contrasts? The basic unit of instantaneous convective rainfall is the convective cell O(10 km), often embedded in meso-scale convective clusters, tropical depressions and/or cyclones 1,2,35-37 . These are sometimes accompanied or followed by more extensive (a) Skill of week3+4 rainfall 24     V. Moron and A.W. Robertson stratiform rains. 38,39 A useful diagnostic of biweekly rainfall scales may thus be the size of daily wet patches (WP), given by the contiguous area receiving at least N mm of rain. Daily aggregates of rainfall will increase the patch-size somewhat through the integration and movement of rainfall within a day, and a threshold of N = 10 mm (WP10) was chosen empirically from the daily rainfall data (Supplementary Fig. 1). A larger daily WP10 may be expected to lead to larger spatial scales of aggregated rainfall amounts and vice versa. A second factor relevant to biweekly rainfall anomaly spatial scale and predictability is the mean rainfall intensity. Since heavy convective rainfall events tend to be highly localized, a large mean intensity associated with such events will tend to be associated with smaller spatial scales of S2S rainfall anomalies, and potentially lower prediction skill. 31,40,41 Finally, the essential source of S2S rainfall predictability is the impact of S2S phenomena, that modulate or "synchronize" small-scale instantaneous rainfall in space and time, such as the MJO 10,42 and other convectively coupled equatorial waves, [43][44][45] as well as the impact of slowly evolving surface boundary conditions including soil moisture and sea surface temperature (SST) anomalies, most prominently driven by ENSO. 12,46 Figure 3 shows the geographical distributions of these interrelated controls on S2S rainfall predictability, estimated from observed data (see Supplementary Table 1 for their PCs). The 90th percentile of daily WP10 shows a large maximum across the western and central equatorial Pacific, with very small values over Africa and the mean subsidence regions in the subtropics and eastern Pacific, linked to descending branches of the Hadley and Walker circulations (Fig. 3a). The spatial pattern of mean intensity (Fig. 3b) is similar to that of mean rainfall amount (not shown) but with an important difference: the mean intensity tends to peak over the landmasses, especially South and southeast Asia, northern Australia and the Maritime Continent as well as most of South and Central America, while Tropical Africa is largely the exception, even if rather large values are observed over Madagascar and Mozambique, partly due to cyclonic events there. The synchronizing impact of S2S phenomena is broken down by timescale into interannual, and the 20-90 day and 7-20 day intraseasonal bands in Fig. 3c-e, in terms of the fraction of rainfall variance in each frequency band. Interannual variations ( Fig. 3c) peak over the core ENSO region, and are generally larger over the oceans; smallest values occur over equatorial and northern Africa and western South America as well as the regions of mean subsidence. The intraseasonal variances (Fig. 3d, e) both exhibit large values over the subtropical oceans, likely associated with equatorial Rossby waves, [43][44][45] while the 20-90 band also peaks over the equatorial Indian Ocean and western tropical Pacific MJO core region. [8][9][10]42 Over land, the intraseasonal bands have moderate to strong impacts on South Asia, northern Australia, eastern and southern Africa, central America and eastern Brazil (Fig. 3d). Lastly, the <7 days variance provides a measure of weather noise which includes the impact of synoptic transient perturbations as well as the fastest CCEWs, 44,45 and clearly shows a sea-land contrast with highest values over the landmasses (Fig. 3f). This is especially clear over Africa, South and SE Asia, the Maritime Continent and Australia. The values peak above 50 % of the total variance over equatorial and northern equatorial Africa and central and western South America.
A geographical classification of rainfall predictability characteristics We next classify each land or ocean gridpoint into six geographical "S2S predictability types" according to the values  of its seven observed rainfall characteristics mapped in Figs. 2, 3. This classification was done using a fuzzy k-means cluster 47,48 analysis, carried out for 2100 land and 6387 ocean points separately ( Fig. 4 and Supplementary Fig. 2); no spatial location data was used. The six resulting predictability types (clusters of gridpoints on Fig. 4a and Supplementary Fig. 2), derived only using observed rainfall characteristics, are ordered by the ECMWF model forecast skill at those gridpoints (Fig. 1a) which was not used in the cluster analysis. It is striking that a classification of locations based on observed quantities only, without any explicit regression on various sources of skill (e.g., MJO time sequence or ENSO states) leads to a relevant geography of the model's skill (Fig. 1a). Note also that the clustering is rather robust versus the choice the included variables ( Supplementary Figs. 3, 4). Removing either the fraction of interannual variance, or that of the fraction of subweekly variance, or considering only the three bands of interannual and both intra-seasonal bandwidths lead to a similar geography of the S2S predictability types. Highest predictability over land (Type 6; red/pink) occurs over northern Australia, the Maritime Continent, areas of eastern and NE Brazil and few grid-points over coastal Ecuador and Galapagos. Interannual and 20-90 day ISO synchronization are largest for this type, highlighting the role of ENSO and the MJO, and spatial scales and daily wet-patch areas are both very large. Perhaps surprisingly, mean daily intensity is above the continental average, (horizontal dashed line), but this is counteracted by the larger spatial scales and S2S synchronization.
Land areas with moderate S2S skill are split into three clusters (Types 3-5). These extend over large parts of southern and eastern Africa, South and Central America, and subtropical Australia. The skill in Type 3 may be increased through relatively large WP10 while Types 4 and 5 exhibit relatively large synchronization at interannual and intraseasonal time scales.
Predictability types 1 and 2 have the lowest skill, located over large parts of west and central Africa, except coastal areas (Type 1), and most of South and SE Asia (except NW and southern India), SE Brazil and Uruguay and parts of central America, SE Africa and Madagascar (Type 2) These low predictability types are quite distinct: Type 2 has high ISO synchronization (especially 7-20 day), as well as slightly larger WP10 than the land average, all of which may increase predictability-but the mean intensity is very large, tending to decrease it. On the other hand, Type 1 (over Africa) lacks interannual or intraseasonal synchronization, and rainfall scales are very small. The impact of <7-day variance clearly peaks for this cluster.
The ocean rainfall predictability types ( Supplementary Fig. 2) are organized around the central Pacific, where the skill is the highest, in almost concentric areas. Note the strong impact of 20-90-day band on cluster 3. The skill is lower across the Atlantic basin, but also over subtropical northern Pacific and western and southern Indian Ocean (Types 1 and 2), mostly in association with small WP10 and a poor synchronization at interannual time scale ( Supplementary Fig. 2).

DISCUSSION
The empirically derived S2S rainfall predictability types over land demonstrate a one-to-one agreement between the spatial scales of observed biweekly rainfall anomalies and ECMWF model week 3 + 4 forecast skill. This is a remarkable agreement of empirical and model results. While models can be expected to improve in the future, this result underlines the expectation that S2S rainfall skill is fundamentally limited in several regions of the tropics, which either lack the predictable climate signals in the S2S band, or in which the inherent scales of rainfall are too small, or rainfall intensities are too large, to make biweekly averages usefully predictable. Most prominently, large parts of West and Central Africa and the Sahel (Type 1) lack S2S rainfall predictability, on average, due to the absence of strong intraseasonal and interannual dynamical controls, coupled with the small spatial scales of daily rainfall; the latter is consistent with the short duration of the convective (or MCS) events across this area 49 and the large amplitude of the diurnal cycle. 50 Note that the decadal time scale, which may be related with significant SST forcing across West Africa 51 is not sampled enough by the period analyzed here from 1996.
In contrast, the other broad region of low S2S rainfall predictability (Type 2) over South and SE Asia does experience temporal synchronization of rainfall from both intraseasonal bands, but the mean intensity (and thus noise) is very large. 31,40 Note also that WP size and mean rainfall intensity are positively related (viz. Figs. 4b, c), consistent with the scaling properties of tropical rainfall, at least over the ocean. 52 The impact of WP10 area on the signal-to-noise ratio thus represents a trade-off between the spatial integration (increasing the signal) and the maximum rainfall received near the center (increasing the noise). Improved and higher resolution models may improve skill here since dynamical sources of predictability and skill are present, as may the consideration of rainfall frequency as the predictand variable. The role of the mean intensity suggests that S2S predictability may decrease due to global warming where a general increase of rainfall intensity in wet areas/seasons is expected. 53,54 Regions of higher continental skill over SE Asia and Australia (Type 6), eastern and southern Africa (Types 4 and 5) and South America (Types 3 and 6) are mostly related to an increased amplitude of the temporal synchronization due to inteaseasonal and interannual bands, and where the relative influence of synoptic time scales <7 days is strongly reduced (over southern Africa). The locations of Types 5 and 6 across the landmasses (Fig.  4a) coincide rather well with areas where either seasonal predictions or impact of MJO have been established. 9,13,23,29,55,56 One caveat to the above classification is that the situation over the continents is complex and that the contribution of the analyzed drivers is only partial, as suggested by the intra-cluster spread in Fig. 4j. Our results pertain to the entire calendar year and predictability may be higher in specific parts of the calendar year, and especially during the early and late phases of the monsoon seasons including the onset and/or demise dates. 57,58 Over the oceans, the skill is larger and the contribution of the S2S climate drivers is stronger than over land. The predictability types are organized broadly in a symmetric way, around the tropical Indian and Pacific Ocean where the skill/scale peak. The large and strong dynamical control of interannual (central Pacific) and mostly 20-90-day (western Pacific and Indian ocean) increase the skill/scales there. As for Type 1 over the landmasses, Types 1 and 2 over the oceans (well extended over the tropical Atlantic and subtropical margins elsewhere) are penalized by the small size of the wet patches while the synchronization provided by interannual and 20-90-days is weak to neutral.

DATA AND METHODS
Daily rainfall data Rainfall amount was extracted from the GPCP-1DD product (resolution 1°× 1°, daily, version 1.3) from 1 October 1996 to 31 December 2017. The GPCP dataset blends various satellite estimates and the daily amounts are calibrated with raingauges. 59,60 The horizontal resolution already filters some of the noise since the minimal size of wet patches or observed spatial scales is roughly 11000 km 2 , which is far larger than the basic unit of convection (~10-100 km 2 ) and corresponds broadly to the lower bound of the meso-scale convective systems (MCS). Thus, we are unable to say if a wet cell is related to a single localized thunderstorm or an amalgamation provided by an MCS. The daily resolution is a limit of the analysis, but wet events usually last for more than one day. 2,49,50 The calibration using rain-gauges, mostly located across the landmasses, may be a source of heterogeneity in the spatial variations of the characteristics defined using observed rainfall. But spatial variations of either WP10, spatial scales of biweekly rainfall anomalies, and fraction of bandwidth variations are very similar in uncalibrated Outgoing Longwave Radiation with the same resolution (not shown).
Daily ECMWF rainfall Daily ECMWF model rainfall was obtained from an 11-member ensemble of S2S reforecasts for the 1998-2017 period (105 starting V. Moron and A.W. Robertson dates each year). The ECMWF forecasting system is a merged 47-day ensemble system updated two times a week that has been run since 2002. 15 The ECMWF runs are performed "on the fly", with a 20-year reforecast set of runs made every time that a real-time forecast is issued. We use the model configuration, which was operational in 2018 (version CY43R3 for the forecast-issue dates from 1 January to 6 June and version CY45R1 for the starting dates from June 6 to December 31). The differences between both configurations are minor and are assumed to be negligible for the analyses performed in this study. Both versions use an atmospheric model at Tco639 (about 16 km) up to day 15 and Tco319 (about 32 km) horizontal resolution after day 15 with 91 vertical levels. The atmospheric model includes land surface interactions and is fully coupled to ocean and sea-ice models. The sea surface temperatures and continental surface characteristics are realistic at the starting dates. Daily rainfall accumulations are available on a 1.5°grid and the week 3 + 4 rainfall is computed from days 15 to 28. The skill is computed considering the same dates between the ensemble mean and GPCP-1DD rainfall after having linearly interpolated the GPCP-1DD rainfall onto the 1.5°grid. The potential predictability (PP) is defined as the mean Spearman correlation between anomalies of each run vs the ensemble mean. As for the GPCP-1DD rainfall, ranks are used instead of amount themselves and slots when the climatological mean is lower than 1 mm/day (from the climatology of the ensemble) are excluded from the computations.
Defining the spatial scales and wet patches areas The daily rainfall data from GPCP-1DD are used to compute the spatial scales of 15-day rainfall anomalies (Fig. 2) and various characteristics of the rainfall variability (Fig. 3). The spatial scales are computed as the area where the correlations are above 1/e (~0.3679) around each grid-point. 49 The rainfall anomalies are computed as the differences between the running 15-day amounts and the climatological 15-day average. The periods when the climatological daily average (daily mean smoothed by a 1/90 c-p-d low-pass recursive filter) are not considered and the anomalies are converted into ranks to cancel the effect of very large positive anomalies. So, the spatial scales are estimated using Spearman's instead of Pearson's correlations. The wet patches are also defined from the GPCP-1DD fields. Wet patches (WP) are defined as the contiguous area (either by gridbox edge or corner) receiving at least 10 mm (WP10). A buffer zone of 10°is used to define the WPs on the 30°N-30°S domain. Considering a lower threshold of 1 or 5 mm leads to patchy, but very large, structures while higher thresholds usually limit the size of most of WP to few grid-points. A higher threshold tends to give very small samples for some areas. There are a total of~730000 WP10 and more than 30% cover only a single grid-point. To compute the 90th percentiles of WP10, each WP10 area is attributed to each of the grid-points belonging to it, and the percentile is computed locally. In that context, a given WP10 is counted more than once as soon as it covers at least two contiguous grid-points.

Daily mean intensity
The mean daily intensity is simply computed as the total amount of rainfall across the whole set of days divided by the total frequency of wet days ≥1 mm. This characteristics mixes the instantaneous rain rate, that may depend mostly on the intrinsic nature (i.e. convective vs stratiform) of the rain-bearing clouds, and the intra-day duration of the wet event that may be mostly driven by the dynamical control of the rainy system (a longer duration is expected for MCS-TD system vs purely local thunderstorm).
Temporal decomposition of the rainfall variance The fraction of variance belonging to four different frequency bands (Fig. 3) is computed from square-rooted daily rainfall.
Anomalies are computed from the climatological daily mean from these square-rooted daily rainfalls, and a recursive filter is used to extract the low-pass (>365 day), band-pass (20-90 day and 7-20 day) and high-pass (<7 days) variations. The 4 bands are expressed as their variance divided by the unfiltered variance. The sum of the 4 bands is usually close to 1.
Clustering of the rainfall drivers of the skill Both a hard and fuzzy k-means cluster analysis 47,48 are applied to the standardized anomalies of the 7 observed variables shown in Figs. 2, 3, separately for the 2100 (land) and 6387 (sea) grid-points. Dry areas which never reach a climatological daily mean of 1 mm/day are not considered in the computations. A pre-processing is done with a principal components analysis, retaining 90% of the variance in the four leading PCs. The goal of this clustering is to give a sound and interpretable regionalization (Fig. 4a and Supplementary Fig. 2a), explaining at least some of the spatial variation of the skill of Fig. 1a,  4j and Supplementary Fig. 2j). The hard k-means algorithm partitions the data into k clusters by iteratively maximizing the ratio of the sum of the between-group variance vs total variance, with the entire process repeated 1000 times from random initial seeds. The 2100 land grid points (and 6387 over ocean) are divided into k distinct clusters, where each gridpoint point can only belong to one cluster. A relatively large increase of this ratio occurs for k = 3 (land), k = 6 (land and ocean), and then for k = 13 (land) and k = 9 (ocean) (not shown). The latter solutions leads to far more complicated maps (not shown) than Fig. 4a and Supplementary Fig. 2a. A fuzzy k-means with k = 6 is then applied to the four principal components explaining 90% of the variance for land and ocean separately. When the value of the fuzzifier is >2, all memberships tend to 1/6, and thus a meaningless clustering, while when it is <1.5, it tends progressively to hard clusters, with all memberships tending to 1. The subjectively chosen value of 1.75 is thus relevant to distinguish core and marginal areas belonging to the six clusters.
The 7 variables included in the clustering are not independent to each other, e.g. the fraction of variability conveyed by the four complementary frequency bands are related by construction. However the pattern correlations between them (Supplementary Table 1) are mostly moderate, suggesting that they may thus contribute differently to the building of scales and skill of timeaveraged rainfall anomalies. Nevertheless, the sensitivity of the clustering to the included variables should be tested. We present three alternative clusterings in Supplementary Fig. 3b-d for land and in Supplementary Fig. 4b Figs. 3b, 4b), 92% of grid-points over land and sea belong to the same cluster; without high frequency variations, 90% (79%) of grid-points over land (sea) (Supplementary Figs. 3c, 4c) belong to the same cluster. The order of the clusters from the point of view of spatially averaged skill used to rank the six types on Fig. 4a and Supplementary Fig. 2a sometimes switches between contiguous clusters (for example between clusters #1 and #2 or between clusters #3 and #4). The third sensitivity test is more stringent by considering only the 3 slower bands of variations in the clustering (Supplementary Figs.  3d, 4d). Even in this case, where the explicit measures of noise (high frequency <7 days and mean intensity) are excluded, the geographical pattern as well as the relative skill rankings (revealed by the colors from blue/purple for low skill to orange/red for high skill), are broadly similar. In fact, 58% (land) and 55% (sea) of the grid-points are in the same clusters as Fig. 4a and Supplementary  Fig. 2a. Nevertheless, the map is slightly noisier and small pockets of different clusters within the regionally consistent ones appears, V. Moron and A.W. Robertson such as over SE Asia and Central Amazonia. Overall, the sensitivity analysis results confirm that the clustering is rather robust and support including the seven variables as different facets of noise and signal (Supplementary Table 1). It is also possible that other relevant variables may be included, but is beyond the scope of this study.

CODE AVAILABILITY
The statistical analyses have been made using Matlab software and the corresponding codes that support the findings of this study are available upon request from the corresponding author.