Extreme climatic characteristics near the coastline of the southeast region of Brazil in the last 40 years

The southeastern Brazilian coast is a vulnerable region to the development of severe storms, mainly caused by the passage of cold fronts and extratropical cyclones. In the last decades, there has been an increase in the occurrence of subtropical cyclones. This study investigates trends and climatic variations, analyzing surface meteoceanographic series at six grid points from the reanalysis databases of ERA-Interim and ERA5 (European Center for Medium-Range Weather Forecasts (ECMWF)) from 1979 to 2018 over the ocean region bounded, approximately, at 18°S, 25°S, and 37ºW, 45ºW (between the states of Espírito Santo, Rio de Janeiro, and São Paulo). Non-parametric statistical tests and the generalized extreme value distribution are employed for annual, seasonal, and daily maxima/minima. The numbers of occurrence of extreme values, as well as the extremal index, are also estimated in order to better understand the behavior of extremes. Annual maximum sea-surface temperature anomalies of the ERA-Interim databases show very low negative values, mainly at the beginning of measurements (between 1979 and 1982), leading to high positive trend values. The results are compared to the updated data from ERA5 which have anomalies that are more homogeneous with positive trends but without statistical significance. The other meteorological series of the ERA-Interim does not present discrepancies. Only the maximum anomalies of air temperature have significant annual and seasonal positive trends at grid points near the coast of Rio de Janeiro and São Paulo. Despite that the analyses for pressure and wind speed anomalies do not indicate significant trends, they present increases in the interdecadal pattern of the numbers of occurrence of extreme percentiles for almost every grid point. Return levels for 10, 25, 50, 75, and 100 years are estimated at each grid point, and many maximum/minimum peaks are close to the return levels for 100-year return periods. The extremal index suggests average cluster sizes associated with no predominance of clustering for the extreme percentiles, which represents weak dependence between the exceedances. These results characterize some independence between extreme meteorological events such as the event that has been taking place in the region. The occurrence of maximum daily wind speed peaks calculated in austral spring, whose values exceeded the previous ones, is identified at three grid points near the southeast Brazilian coast, caused by the passage of the subtropical cyclone “Deni,” which occurred in November 2016.


Introduction
The study of extreme climatic patterns has become increasingly relevant because of the increase of frequency of severe meteorological systems that cause destructive impacts when combined with several factors such as disorderly growth of cities and increase of urban population (Easterling et al. 2007). Air pollutants and marine contaminants contribute to the environmental changes (IPCC 2009). The number of natural disasters that have been occurring since the middle of the twentieth century and the beginning of the twenty-first century is mainly attributed to the increase in the concentration of greenhouse gases, and according to Ryan (2015) and Di Giulio et al. (2019), it is necessary to develop governance strategies to mitigate their effects. Wallemacq (2018) shows that between 1998 and 2017, thousands of people have suffered from extreme weather events with 91% of all disasters being related to droughts, floods, and storms. Therefore, climate change impact in recent decades is an important topic to be analyzed for improving our understanding of the causes and consequences with regard to global warming (D'Agostino and Lionello 2017; Tegegne et al. 2020).
An important aspect according to Collins et al. (2018) is the development of meteoceanographic systems resulting from the warming of ocean and atmosphere. This development is related to the increase of evaporation and humidity, leading to the occurrence of extreme events (Xie 2020;Pereira et al. 2020). The heat balance in the ocean-atmosphere interaction can lead to small variations in sea-surface temperature (SST), which can lead to large variations in heat flows, playing a relevant role in local, regional, and global climatic conditions (Pezzi et al. 2016). Cheng et al. (2019) found that the excess heat stored in the oceans has been causing temperature increases. Several studies in recent decades have shown positive trends in SST, which may result in changes in the pattern of ocean-atmosphere interaction (Bombardi et al. 2014;Turner et al. 2016). Cheng et al. (2020) showed an SST trend above 1.2 ºC/100 years in the Atlantic Ocean between 50ºS and 30ºS latitude.
Significant meteorological events, in South America, have been registered, and the highest rate of cyclone formation, in the Southern Hemisphere, is seen between the Antarctic Peninsula and the Southern region of Brazil (Hoskins and Hodjes 2005). Transient atmospheric phenomena on a synoptic scale, described by these authors as being related to oceanic variability in the South and Southeast regions of Brazil, account for some cyclogenesis mechanisms. Previous studies on precipitation variability and trends in SST in the southeastern region of Brazil have shown an increase in SST in recent decades (Silva Dias et al. 2013).
The southern and southeastern Brazilian coasts are vulnerable regions where the development of storms in the ocean area occurs. The most frequent cyclones are extratropical ones that develop over well-defined oceanic regions of southeastern South America, moving with cold fronts (Gan and Rao 1991;Reboita et al. 2010;de Oliveira et al. 2011). However, in the last decades, an increase in the frequency of occurrence of subtropical cyclones has been observed near the coastline of the southeastern region of Brazil (Gozzo et al. 2017;Reboita et al. 2019). According to Evans and Braun (2012), 1.2 subtropical cyclones per year occur, on average, over the South Atlantic basin with little monthly variability. Gozzo et al. (2014) simulated climatology of subtropical cyclones and their synoptic pattern over the South Atlantic, proposing a broader definition of these systems. The authors concluded that most of the subtropical cyclones develop during the austral summer (December to February), but in autumn (March to May), the tropical environment leads to the most intense episodes resulting from seasonal meteorological conditions. Because of this increase in the frequency of subtropical cyclones in the Brazilian maritime region, regulations for marine meteorology activities were approved to monitor and classify these systems through NORMAM 19-2011 (Maritime Authority Standards for Maritime Meteorological Activities, Brazilian Navy). This normalization was accessed from https:// www. marin ha. mil. br/ dhn/?q= en/ node/ 264 and follows the terminology adopted in the World Meteorological Organization (WMO), No. 1194/2017, Global Guide to Tropical Cyclone Forecasting accessed from https:// cyclo ne. wmo. int/ pdf/ Global-Guide-to-Tropi cal-Cyclo ne-Forec asting. pdf.
Therefore, the interest in this research is to better understand the behavior of surface meteoceanographic variables over the ocean area near the coastline of the southeast region, which has been affected by these extreme events with more frequency than previously recorded, causing material and human damage. Their occurrence is more common in the southern region of the country. However, the study region is an area with limited temporal and spatial data available for long-term analysis. Thus, we use reanalysis data (Compo et al. 2011) over the maritime area, covering a spatially extensive region near the coast. Initially, the study is based on choosing grid points for the analysis of extreme values over the southeastern coast of Brazil. After several searches, the extreme events at some grid points during the studied period are identified, making it possible to obtain new information about their occurrence. A methodology focused on statistical analysis is used to describe the temporal variations as well as on the analysis of extremes for a long period.
The purpose of this study is to first apply non-parametric statistical analyses of Mann (1945), Kendall (1975), Sen (1968), andPettitt (1979) to test stationarity and eventual abrupt changes in time series. The magnitude of relative changes and significant trends are identified during the last 40 years. In addition, extreme value analysis (EVA) is also applied in order to better understand the behavior of these variables, as well as the probability of occurrence of extreme events in the region. This approach has several applications for environmental data (Gilleland and Katz 2006;Cheng et al. 2014;Mascaro 2018). The generalized extreme value (GEV) family of distributions is used because of the flexibility in fitting the particular cases of the distributions relative to the maximum or minimum blocks of the time series according to the definitions described by Fisher and Tippet (1928). The extremal index is also used to measure the degree of dependence/independence of several extreme observation sequences, modeling clusters of extremes within a short period of time (Coles 2007) and so to know better the behavior of extremes during the studied period.

Study area
The study area is located over the oceanic area of the southeast region of Brazil between 18ºS and 25ºS and 37ºW and 46ºW off the coast. A gridded data set was chosen with a spatial resolution of about 778 and 875 km between the latitudes and longitudes, respectively, making it possible to detect the presence of synoptic scale phenomena (Orlansky 1975). Figure 1 shows the geographical location of the study area with the six selected grid points from the reanalysis ERA-Interim and its update as ERA5 databases over the South Atlantic Ocean basin.
The atmospheric circulation pattern over the study area is related to the general circulation of the South Atlantic under the domain of the South Atlantic subtropical high (SASH) that is a semi-stationary anticyclone system, which presents a seasonal movement. During the summer, this system becomes weaker than winter, moving southerly, and on its western side, the winds blow northeasterly towards the southeastern coast of South America (Fig. 2), being more intense in the southeastern coast of Brazil (de Oliveira et al. 2012). Periodically, this subtropical anticyclone is affected by the passage of frontal systems that move from the southwest across the northeast along the Brazilian coast (Campos et al. 1995). Therefore, this region is strongly influenced by the wind pattern that influences the distribution of water masses classified into three categories according to temperature and salinity (Piola et al. 2000;Soares and Moller 2001), such as tropical water (TW), coastal water (CW), and South Atlantic Central Water (SACW). Figure 3 show the dynamic along the east coast of Brazil is characterized by the presence of the South Equatorial Current that flows towards the south and forms on the west side of the subtropical gyre of the South Atlantic called the Brazil Current (BC), originating from the bifurcation of the Equatorial Current, around 10ºS. The western border of the subtropical gyre is parallel to the South American continent, characterizing a medium circulation that is distributed in a closed anticyclonic gyre (Stramma and England 1999;Piola et al. 2000). The BC flows southward along the continental shelf, reaching the Argentinian coast at approximately 38ºS where the Brazil-Malvinas Confluence (BMC) is located. Thereby, the convergence between the BC and Malvinas Current (MC) results in the generation of meander and intense vorticity (Peterson and Stramma 1991), leading to the development of intense cyclogenetic activities, such as extratropical cyclones. The action of the winds on the sea surface has spatial similarities with the subtropical atmospheric gyre over the South Atlantic basin (Garzoli and Matano 2011), called the South Atlantic subtropical high (SASH), which plays an important role in the regional climate over southeast South America (Bombardi et al. 2014;Sun et al. 2017).

Meteoceanographic databases on the regional scale
ERA-Interim reanalysis data of the European Center for Medium-Range Weather Forecasts (ECMWF) are used for the 40-year-long period from January 1979 to December 2018 at 6 hourly intervals (00, 06, 12, and 18 UTC) of mean sea-level pressure (MSLP), zonal (u), and meridional (v) wind components at 10 m height, 2-m air temperature (T), and sea-surface temperature (SST). Six grid points were chosen along the continental shelf of the southeastern Brazilian coast between 18ºS and 25ºS and 37ºW and 46ºW off the coast (Fig. 1). The wind speeds (WS) were calculated from the wind components (u and v) for each selected grid point. However, during the statistical analyses of SST, it became necessary to compare these data with other data at the same points and period, because of discrepancies found at the beginning of the observations. This comparison was a difficult task because the study region is comprised in the oceanic area with a sparse data network. Thus, it was decided to use the same source of the ECMWF whose SST data had already been updated as ERA5 since 2019, downloaded for the same period and same grid points.
Originally, the ERA-Interim started in 1989, but the 10-year extension from 1979 to 1988 was produced in 2011 and continued in real time until 2019, after which an update was performed as ERA5 (Hersbach et al. 2020) with a grid of 1.0º × 1.0º latitude/longitude. Comprehensive documentation of the ERA-Interim system can be found in Berrisford et al. (2011) and Dee et al. (2011). Further information on the ECMWF data assimilation system can be found in the online documentation available from http:// www. ecmwf. int.
According to Hirahara et al. (2016), the SST data set from ERA-Interim is a combination of five products that generate a single series from 1979 to 2019 (Reynolds et al. 2002;Reynolds and Chelton 2010). These authors evaluated the latest generation of high-resolution analyses for SSTs through intercomparison between products and through the assimilation of atmospheric data, leading a priori to eliminate uncertainties in the reanalysis and proposed the use of SST through ERA5 (Donlon et al. 2012;Hersbach et al. 2020). The analyses, in most oceanic regions, use interpolations and adapt satellite data to data in situ (Reynolds et al. 2002) because of sparse data in these regions. The advanced veryhigh-resolution radiometer (AVHRR) has greater coverage and longer registration. According to McClain et al. (1985), the SST has been estimated with the aid of these instruments since 1981. The assimilation of radiance data, which is affected by systematic errors in the radioactive transfer models, is embedded in the system (Dee and Uppala 2009). Therefore, AVHRR SSTs are, in principle, equivalent to measurements in situ (although measurements away from the buoys used in the regression cause some degree of uncertainty); however, they cannot be used in combination with in situ data without adjustment (Reynolds and Chelton 2010).

Mann-Kendall trend test, Sen's slope estimator, and Pettitt's test
In this research, non-parametric methods are used to detect trends, which make no assumptions about the distribution that the data follow, requiring only independence of the sample data. Several methods have been developed to detect changes or discontinuity on climatic series over time from standard statistical techniques. The most commonly used techniques apply Student's or Mann-Kendall tests with a certain degree of confidence, requiring enough data of around 10 or more years (Rodionov 2004).
The application of the Mann-Kendall (MK) test is recommended by the World Meteorological Organization (WMO) to detect statistically significant long-term monotonic trends in environmental data (Gavrilov et al. 2016;Kocsis et al. 2020). The non-parametric MK and Pettitt tests are the most used to detect both the lack of homogeneity and the changepoint in the mean of a time series (Arikan and Kahya 2019).
The tests are used here on maximum and minimum annual values, which are therefore considered independent values over time. Seasonality and trend are not excluded from the data because of the statistical analyses applied to search for any temporal changes, and it is generally contraindicated to remove trends/seasonality before doing EVA but rather to model those elements.
The null hypothesis ( H 0 ) of the MK method is that there is no trend in the time series and the alternative hypothesis ( H a ) is that there is a trend for a given level of significance p. H 0 is not found to be significant if p is greater than a pre-chosen significance level, , chosen here to be 5%. The result of the analysis is the Tau score, indicating a decrease or increase in the attributes of the series when it is negative or positive, respectively (Appiotti et al. 2014). Although efficient and robust, the MK test does not provide trend values, but it can be complemented by the non-parametric slope estimator proposed by Sen (1968). This method fits the slope to the sampling points by choosing the median of the slopes, considering pairs of points that vary per unit of time.
It uses a linear model to estimate the slope and variation of the residuals that must be constant over time (Tao et al. 2014;da Silva et al. 2015). However, to detect change in the mean of the series with statistical significance (Back 2001;Liu et al. 2012), a version of the Mann-Whitney test uses Pettitt's test where H 0 considers homogeneity in the series; otherwise, H a suggests that a change-point or breakpoint occurs. These tests are applied to annual maximum T, SST, and WS and annual minimum MSLP, as well as to seasonal maximum (minimum) values, using the 5% statistical significance level.
As the interest of this research is focused on the behavior of extreme surface variables to better understand any changes in climatic patterns, we investigate the number of occurrence (Noc) above (below) thresholds considered extreme percentiles chosen according to results found by Reboita et al. (2019). These results refer to the range of each surface variable in the development of subtropical cyclones. These authors documented the main features of six subtropical cyclones that occurred between 2010 and 2016 over the southwestern South Atlantic Ocean near the Brazilian coast. Most subtropical cyclones present MSLP ranging from 1005 to 1010 hPa and maximum 10-m WS ranging between 10 and 15 m/s in the initial phase. This study also showed these systems develop over a wide range of SST (from 22.5 to 28.6 ºC), indicating a weak dependence of the cyclogenesis on the warm SST (Guishard 2009;Gozzo et al. 2014) instead of the values of thresholds (26 ºC to 26.5 ºC) normally favorable to tropical cyclones (Trenberth and Shea 2006).

Application of extreme value analysis
The basic elements of extreme value theory (EVT) were initially developed by Fisher and Tippett (1928) and Gumbel (1954) to investigate events that are considered extreme or rare, enabling the behavior of these events to be described and quantified by estimating a possible limit distribution for the maximum (minimum) of the sample composed of independent and identically distributed (iid) random variables. These authors proved that the three possible types of asymptotic distribution of extreme values, known as the Gumbel, Fréchet, and Weibull, arise whenever such a limiting distribution is non-degenerate. In other words, if a limiting distribution exists for the maxima of a series, then it must be one of these distributions. The behavior of the tail of F(x) determines this condition with a large data set, and, therefore, the asymptotic properties of the extremes can be determined without knowing the properties of the random variables in detail (Leadbetter et al. 1983;Coles 2007). Extreme value analysis (EVA) allows for the evaluation of the probability of occurrence of more extreme events from sample sequences than previously observed (Gilleland 2015;Towler et al. 2020). The theory states that under certain regularity conditions, it is possible to estimate the likelihood of the occurrence of extreme values based on a few basic assumptions. EVA refers to analyzing events with a low probability of occurrence.
The random variables taken over suitably large (e.g., annual) block maxima, which have a non-degenerate distribution, follow the generalized extreme value (GEV) distribution function (df) that encompasses the three types mentioned above (Perrin et al. 2006;Coles 2007). Similar theory gives that the generalized Pareto distribution (GPD) is the limiting df for excesses above a high threshold provided the limit exists. Both the GEV and GP can be characterized in terms of a Poisson process, which allows the parameters related to the frequency and intensity of extreme events are fitted simultaneously (Gilleland and Katz 2016).The GP and PP approaches, collectively known as the peak-overthreshold (POT) approach, have the advantage of using more data but often at the expense of violating the assumption of independence between excesses. If there is enough data to fit a GEV to the block maximum/minimum, it can be a cleaner and simpler approach, as the assumption of independence is more likely to be met and its quantiles exactly correspond to the return levels. Therefore, to estimate the parameters of the GEV df (Gumbel, Weibull, and Frechét), it is necessary to fit it to a set of maximum or minimum values of a long data series (≥ 30 years) recorded each year (Malevergne et al. 2006;Escalante-Sandoval 2013).
The analysis of extreme climatic data is important for better understanding severe atmospheric processes for helping the evaluation of possible environmental risk (Tegegne et al. 2020). One of the important characteristics in studies of extreme environmental values refers to the distribution and tail size of the distribution where the limit value is located (Mascaro 2018;Telesca et al. 2020). The three types of distributions encompassed by the GEV family are associated with the value of a shape parameter (ξ) (Coles 2007). The GEV family of distributions is shown in Eqs. 1 and 2: where: σ > 0 (scale parameter). ξ ∈ ℟ (shape parameter). Type I: Gumbel (light tail, ξ = 0) domain of attraction for many common distributions.
Type II: Fréchet (heavy tail, ξ > 0). Type III: Weibull (bounded upper tail, ξ < 0). There are different methods available for performing parameter estimation of the distribution (location, scale, and shape). The maximum likelihood estimation (MLE) method is used here because of its wide applicability and efficiency (Escalante-Sandoval 2013) to evaluate the suitability of the GEV hypothesis in estimating the parameters and modeling the annual extremes. This method allows for estimating the GEV parameters through hypothesis tests related to linear and non-linear restrictions to the parameter vector. Its theoretical justification provides more security to the methods used in risk management (Escalante-Sandoval 2013; Benstock and Cegla 2017). Hence, for samples composed of independent random variables of a univariate density function, the likelihood function (Coles 2007) is where Ψ is the parameter vector to be estimated and f (x i ;Ψ) is the univariate probability distribution function (pdf).
The logarithm of the likelihood function is often more convenient to solve for , and thus Eq. (3) becomes, as the log-likelihood or when ξ ≠ 0 Analyzing changes in the frequency or intensity of extremes, which are further away from the tail, is associated with more uncertain results because of their low frequency of occurrence. Percentile-based extremes are often used in the analysis of climatic variables because of the increase in the sample size that allows for better statistical assessments. There is great interest in estimating the extreme quantiles through the GEV distribution because of their interpretation as return levels associated with the 1/p return period, where p is usually given in years (Benstock and Cegla 2017;Telesca et al. 2020). Return levels are calculated according to the assumption that the maximum value is exceeded, on average, once every p years with probability Pr {Y m > Q n , p } = 1/p. The GEV df fit to block maxima of the data where Q n,p is the quantile 1-1/p of the limit GEV df. We seek z p such that G(z p ) = 1 − p, where G is as in Eq. 2, and considering y p = − 1/ ln(1 − p), the associated return level z p is given by Eq. 6: The generalized Pareto (GP) df is a model for the excesses over a high threshold. It is given by where: u is a high threshold (x > u). σ u is the scale parameter σ u > 0. ξ is the shape parameter − ∞ < ξ < ∞. This parameter determines three types of df's: ξ > 0 Pareto (heavy tail); ξ < 0 Beta (upper bound); and exponential as ξ → 0, resulting in Eq. 8: For iid excesses, the return level y R , P is approximately (1 − 1/p) 1/n where n is the number of the quantiles of the excess distribution. However, GP df quantiles cannot be directly interpreted as return levels because the data no longer derives from blocks of equal length. Instead, an estimate of the probability of exceeding the threshold Ҫ u is required. Therefore, the x m value that is exceeded on average once every m observations (i.e., an estimated level of return) is given by Eq. 9: Another statistical tool is the extremal index (ϴ) that is used to identify clusters related to extreme events as well as the degree of dependence between them. The mean cluster size is the number of extreme observations in a cluster where the ϴ value varies between 0 and 1 and can be considered as the fraction of independent events in a selected sample (Brabson and Palutikof 2000). Smaller values of ϴ indicate strong extreme dependence, while ϴ closer to 1 indicates that extremes are independent. Clusters of extreme values result from the dependence in the time series, referring to the occurrence of several extreme records in a short period of time. ϴ is, therefore, a parameter used for identifying such dependence in which successive days of extreme events tend to increase the risks because of the duration of such severe events, measuring the degree of clustering among exceedances (Smith and Weissman 1994;Ferreira 2018;Cai 2019).
To determine clusters of exceedances, it is necessary to choose a threshold u, which defines the exceedances above (6) z p = + y p , for ≠ 0 + lny p for = 0 it. The run length (r) can be used to separate the clusters with consecutive values above a given u, so values below u will not belong to the same cluster, only the next ones above u and so on. A sequence of r extremes can be viewed as a point process in time, indicating the moment at which the extreme event occurs (Coles 2007). A cluster begins once the series exceeds u. The run length, r, is the number of times the series must fall below u consecutively before the end of the cluster is defined and the next time the series exceeds u, a new cluster begins (Coles 2007).
Here, we take an empirical approach to examine different clustering in the sample data for the maximum daily WS (spring) because of the occurrence of extreme WS in the region with the threshold fixed between extreme percentiles. The extremal indices are calculated at each grid point.
The R software package, extRemes, is used to fit the extreme value models to the data (Gilleland and Katz 2006;Gilleland et al. 2013;Gilleland and Katz 2016).

Sea-surface temperature
The differences found between annual maximum SST anomalies of the ERA-Interim and ERA5 databases used in the trend analyses are discussed here, the results of which show discrepancies mainly between 1979 and 1982. Annual maximum SST anomalies of the ERA-Interim show the largest negative variations in 1979 at all grid points. However, in the subsequent years of negative anomalies, these maxima tended to decrease. Non-homogeneity of the series can be identified on trend lines with accentuated slopes resulting from changes in the mean of the series detected in two periods. Initial analysis demonstrates that there is a need for comparison between databases; thus, the ERA5 is used, which is available from 2019, covering the same period. When compared with the ERA-Interim, lower anomaly values are observed, mainly in 1979, showing that the ERA5 series are more homogeneous. Figure S4 shows the annual maximum SST anomalies from 1979 to 2018 from ERA-Interim (HadISST-www. meteo ffice. gov. br) and ERA5 (HadISST-version 2) databases with their respective linear trends. The variability between them can be seen in the graphs by comparing the differences between the anomalies, mainly in the period from 1979 to 1982. The trend lines show that SST Max from ERA5 presents no considerable trends at all grid points, in contrast to ERA-Interim. It is important to note that the negative SST anomalies observed between 1979 and 1982 with the data from ERA-Interim suggest an underestimated measure on the region. Table S1 presents the Pearson's correlation coefficients between the updated SST Max from the ERA dataset for the study region in 1979-1982 period. Higher values are observed in grid points 3(ES) and 4(RJ) and lower values in grid point 6(SP) with an average value of 0.60. For this reason, SST from ERA5 is used, and ERA-Interim is used for the other variables, as they do not show sharp variations.
The negative SST Max anomalies can be associated with several factors such as the flow of cold water (cold tongue) from the MC to the north and the resurgence processes in the region. Positive anomalies may be related to TW towards the south (Pezzi et al. 2016). Barros et al. (2000) showed that positive (negative) SST anomalies between 20°S and 40°S are associated with the movement of the South Atlantic convergence zone (SACZ) towards the south or north. It is more evident in the summer and is characterized by an elongated axis of clouds and winds converging in the northwestsoutheast direction, reaching the southwest Atlantic Ocean (Chaves and Nobre 2004).

Statistical trend tests of the yearly and seasonal meteorological data
The results of the statistical tests are presented only for T Max from ERA-Interim whose trends are statistically significant (p < 0.05) only for points 4(RJ), 5(RJ), and 6(SP). At these grid points, the Tau scores yield low values, but significant increasing trends are observed with Sen's slope around 0.78 ºC, 1.08 ºC, and 1.19 ºC, considering the 40-year period, respectively. The approach of Pettitt's test presented upward change-points between 1997 and 2008 (Table 1). These change-points are showed in Fig. 4 for the time series of annual T Max between the years that the breakpoints are found in each of the three points.
To verify the seasonal trend pattern in T Max , the averages of the maximum values are calculated for each season, considering the following months: January, February, and March (summer); April, May, and June (autumn); July, August, and September (winter); and October, November, and December (spring). The seasons correspond to the oceanic heating cycle of the austral oceanic region where the peaks of temperatures occur in summer. These time periods are used for oceanographic data sets (Appiotti et al. 2014) with a delay around 10 days in relation to the astronomical period. Thus, the averages of the December data do not consider data from the previous year. Table 2 presents the results of the MK tests applied to the seasonal series where the greatest trends are observed for the austral summer (JFM) and spring (OND) also for grid points 4(RJ), 5(RJ), and 6(SP). The results agree with previous studies carried out in the southwestern Atlantic Ocean near the study area (Gozzo et al. 2014), showing that most of the subtropical cyclones develop during austral summer  (JFM), although in autumn (AMJ), the tropical environment leads to the most intense episodes where stronger surface flows and weaker vertical wind shear is more apparent. Figure 5 shows the seasonal anomalies of annual T Max with the highest trends observed during JFM and OND according to Table 2. The anomalies observed during OND show, at all grid points, that there are two distinct periods: one from 1979 to 2006 with negative and low values (around − 0.5 to − 1.0 ºC) and the other from 2006-2008 onwards with positive values around 1.7 to 2.5 ºC, mainly at grid points 5(RJ) and 6(SP), in which intensifications of annual T Max occur. It is important to note that although the analyses were from grid points, taking into account the spatio-temporal nature of the data set, it is possible to identify changes in the studied period, mainly in seasonal anomalies of T. Nevertheless, WS and SST-ERA5 also present positive trends, but without statistical significance, which for these variables may be a result of the sample size. Another relevant aspect is that more intense meteorological systems, such as subtropical cyclones are observed since 2010.     C. Otherwise, for the minimum values, an amplitude of 5º C between grid points 1(ES) and 6(SP) is observed. This amplitude is mainly related to the differences in latitude between the two grid points and because of the passage of migratory systems over the region, as can also be observed for the minima and maxima of the MSLP. The WS maxima also show homogeneity between the points, except at point 5 (RJ) where it is possible to identify the passage of an extreme meteorological event.

Decadal variability basic statistics
The percentile-based threshold defined as a sequence of contiguous values above (below) a given percentile (Telesca et al. 2020) is used to define two thresholds; namely, the 98th percentile in T (> 27 ºC), WS (> 10 m/s), and 5th percentile in MSLP (< 1010 hPa) at each grid point per decade. These thresholds are based on previous studies (Gozzo et al. 2014 andReboita et al. 2019). A characteristic of a run is its length that represents the period m in which the variable under study is above (below) the selected threshold. Hereafter, a wind extreme is defined as a run with a specific length m. Thus, the amount of r occurrences (Noc) above (below) these extreme percentiles is found, and the results of the interdecadal patterns are shown in Fig. 6. It is interesting to note that the Noc of T varies around 150 to 400 values between 1979 and 2008 except in point 6(SP) where there is an increase from 1989 to 1998, reaching around 800 values in the last decade (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). In other points, an increase can be observed from 1999 to 2008, reaching over 1000 occurrences at point 5(RJ) in the following decade. Besides the T decadal variations, increases in the Noc of WS from 1989-1998 are also found at all points; however, in the decade 1999-2008, there is a decrease at points 1(ES) and 6(SP) with a variation between 450 and 300 values. Grid points 4(RJ) and 5(RJ) reach around 700 occurrences in the last decade. For the MSLP, there is more pronounced variation of the Noc at grid point 1(ES); however, a decrease at all grid points from the decade 1999 to 2008 is observed in accordance with the increase seen in the wind series for the same period.

Evaluation of the GEV distribution
The GEV df is fitted to the annual block maxima of SST from ERA5, T, and WS and block minimum of MSLP from ERA-Interim. MLE-fitted parameter values are presented in Table 4 with the respective standard deviation. The values of ξ are in most cases less than 0, indicating the Weibull distribution, confirming previous studies about environmental extreme data as temperature, wind, and others that tend to be modeled more appropriately with bounded upper tails, suggesting that there are finite values that cannot be exceeded  Gilleland 2015). The upper limit is estimated using the formula µ − σ/ξ and, of course, the values of these parameters, and hence, the upper limit must be estimated from the data, and there is much uncertainty in the estimates. The shape parameter is close to zero for WS at grid points 1(ES) through 5(RJ,) and the Gumbel hypothesis cannot be rejected. Table 5 shows the return levels associated with the 10-, 25-, 50-, 75, and 100-year return periods (estimated from the best fitting distribution along with 95% confidence intervals), maximum/minimum values. Note that the return levels for maximum SST, T, and WS and minimum MSLP gradually increase (decrease) for longer return periods, respectively. Analyzing the results, one would expect that these variables will exceed the estimated values, on average, once per each return period, for example, every 10 years up to 100 years. Table 5 also shows that many maximum/minimum peaks are close to the return levels calculated for 100-year return period.
During data analysis at the grid points 3 (ES), 4 (RJ), and 5 (RJ), we identified extreme WS related to the occurrence of a subtropical cyclone between the states of Espírito Santo and Rio de Janeiro in the austral spring (November 2016). It is important to highlight that the analyzed data in these grid points present seasonal variability of anomalies T Max in OND with a positive trend from 2006 to 2008 onwards at the grid points 4(RJ) and 5(RJ) as shown in Fig. 5. Also, the evolution of the interdecadal Noc (Fig. 6) shows an increase in their values, which exceed the extreme percentiles of T and WS from 1999T and WS from -2008T and WS from to 2009T and WS from -2018 decades in this oceanic region. Therefore, results of the quantile graphics for annual T and WS data at grid points 3(ES), 4(RJ), and 5(RJ) are displayed in Fig. 7. The qq plot makes a straight line, suggesting that the assumptions justifying the use the GEV distribution are reasonable for modeling annual extreme data.

Occurrence of subtropical cyclones identified at selected grid points over the ocean area of the southeast region of Brazil
The statistical analyses allow us to identify significant positive trends in seasonal anomalies of T Max (Fig. 5) and in interdecadal Noc of WS and T (Fig. 6), mostly on the These systems produced severe weather conditions such as strong winds, reaching 17 m/s at 10 m height for a long period and high amounts of rainfall along the southeastern coast of Brazil. However, in the period of 40 years analyzed in this work, some severe meteorological systems related to the passage of cold fronts have been identified because they are the most common systems in the region. According to the analysis of percentiles for maximum daily WS for the OND months during the 40-year period, we verify that the 90th percentile concentrates values around 10 m/s. The 90th percentile was then used as a threshold to examine the excesses of maximum daily values above this reference value for the possible development of subtropical cyclones in the study region. Because of the spatial analysis, it is possible to identify the occurrence of a subtropical cyclone with severe characteristics over the maritime area close to the states of Espírito Santo and Rio de Janeiro for at least three grid points where the most significant extreme values are detected in the seasonal maximum daily WS data set. Figure 8 shows the excesses over this threshold at the six grid points, whose peaks are related to the passage of cold fronts associated with the development of low-pressure systems. In Fig. 8b, c, d, and e, it is possible to verify the occurrence of peaks that exceeded the previous ones between October and November 2016. Figure 8c, d, and e show the maximum peaks of WS on November 15, 2016, at grid points 3(ES), 4(RJ), and 5(RJ) related to the subtropical cyclone named "Deni" that occurred between 15 and 17 November 2016.
This event started as a low-pressure system on the continent towards the ocean near the coast of São Paulo and Rio de Janeiro, reaching the stage of a subtropical storm on November 15, 2016. According to Reboita et al. (2019), the development of this cyclone occurred near the south and southeastern Brazilian coast. Wind gusts are recorded around 20 m/s with central pressure around 998 hPa. Its initial position was, approximately, at 23º 00'S/42º 50'W near the state of Rio de Janeiro (Fig. 9) with a 3-day life cycle. Deni caused more precipitation over the ocean than the continent, generating speeds above 30 m/s in its life cycle, impacting the coast of the state of Rio de Janeiro. The extremal index ϴ is therefore calculated for daily maximum WS in spring months using percentile-based extremes as a threshold, in each case, between the 90th and 99th percentile. As the storm length in the region is, approximately, of 1-3 days, the r values between 1 and 3 should be more appropriate, which coincides with the estimated value of r using the algorithm of the extRemes software package.
The range of the WS thresholds is between 9.6 and 13.7 m/s, and the results of θ suggest a mean value around 0.73, which represents weak dependence (values closer to 1) between the exceedances. The average cluster size is around 1.6 to 1.2 for the 90th to 95th and the 96th to 99th percentiles, respectively, showing there is almost no predominance of clustering. We believe that these results agree with the occurrence of rare events such as subtropical cyclones, and they are in accordance with studies developed by Evans and Braun (2012) where the authors verified the occurrence to be on average about 1.2 events per year over the South Atlantic basin with little monthly variability. The choice of r (separation of clusters) yields values of mostly 1, and only grid point 5(RJ) has values between 1 and 5. Thereby, the number of exceedances does not vary with increasing thresholds when compared to the number of clusters ( Fig. 10a and b). We can also observe that the point 6(SP) shows a similar pattern in both graphs to point 2 (ES), where gusts occur with lower intensities, which may possibly be related to the fact that more intense systems did not pass through the chosen grid point in the analyzed period.

Conclusions
The ocean and the atmosphere interact in ways that cannot be considered separately because its interface is closely related to the SST whose oceanic and atmospheric processes are fundamental as sources of potential energy in the cyclogenetic formation and development.
Significant positive trends are observed in the pattern of annual maximum SST in the ERA-Interim data set at all selected grid points over the oceanic region; however, this trend is smoothed when compared with the ERA5 data for the same points and same length series, especially resulting from the intense negative anomalies observed in the initial period of records from ERA-Interim. Therefore, these trends found in SST samples cannot be conclusive. These data should be better evaluated and compared with other sources because the variability in ocean dynamics caused by the increase in SST can lead to the development of extreme climatic systems in the region, a fact that has been observed in recent years.
The GEV model for the reanalysis data shows satisfactory results and suggests the upper-bounded Weibull distribution represents the distribution of their annual maxima. In the grid points 2(ES) and 5(RJ), the Gumbel distribution cannot be ruled out for the wind data. The lower values of the extremal index for the daily maximum WS, considering percentiles around the 90th, do not indicate excesses in consecutive days, which can characterize extreme wind speeds also related to the passage of transient systems, while the most extreme percentiles may be related to the most severe events that have been occurring in the region.
The methodology developed based on the spatial and temporal statistical analysis (1979 to 2018) of the extremes of the meteoceanographic variables allows for the identification of areas over the region of the states of Espírito Santo and Rio de Janeiro where positive trends are occurring, mainly in wind speed and air temperature series. The analysis of the variability in the series of WS during the last 40 years makes it possible to identify the passage of a subtropical cyclone through the selected points near these two states, which is considered a rare event in the region. The extreme peaks that exceeded the other values in the study period agree with the rare event that occurred in this region. Thus, it becomes more appropriate that the spurious (outliers) values observed in the time series should be analyzed more rigorously, since the occurrence of climatic extremes, in general, has been increasing in intensity and frequency in various regions of the planet.
The results from this research are particularly important with regard to projections of climate change that predict increases in the frequency and severity of extreme events in the next decades. Our studies show the importance of better understanding the mechanism on ocean-atmosphere interaction in the southwest Atlantic Ocean, and for that, it should aim to increase the quantity and quality of databases in the oceanic region near the Brazilian coast to identify future climate changes relevant to climate control of the South American continent.

Author contribution
The first author designed the study, methodology, data, and result analysis and wrote the manuscript. Jorge Luiz Fernandes de Oliveira collected samples, is responsible for data processing, and contributed for result analysis. Pedro José Farias Fernandes is involved with the statistical program in R. Eric Gilleland contributed to the analysis of the extRemes software results and corrected the manuscript. Nelson Francisco Favilla Ebecken contributed to formal analysis and reviewed the manuscript. All authors read and approved the final manuscript.
Funding The first author would like to thank the support provided by the Foundation Carlos Chagas Filho Research Support of the state of Rio de Janeiro (FAPERJ).
Code availability v72i08.R: R replication code.

Declarations
Ethics approval Not applicable.

Competing interests
The authors declare no competing interests. Fig. 10 a Number of excesses over WS thresholds and b clusters using GEV model: a excesses over thresholds from 9.6 to 12 m/s 1(ES); from 10.1 to 12.5 m/s 2(ES); from 10.7 to 12.9 m/s 3(ES); from 11.4 to 13.4 m/s 4(RJ); from 11.7 to 13.7 m/s 5(RJ), and from 10.3 to 12.5 m/s 6(SP). b number of clusters after declustering from 238 to 29 clusters in ascending order of u-values with r between 1 and 5