Abstarct: Abstract Statistical analysis of hydrological variables, which is of great scientific and practical importance in the design and operation of water resources systems, is usually carried out baxsed on the assumption of stationarity. Basically, long-term climate change causes changes in climate variables such as precipitation, temperature, and sea level over time. In this thesis, frequency analysis of extreme values, time series of annual maximum values (AMS), and peaks over threshold (POT) of precipitation and sea level variables were extracted. Then, the fitness of different statistical distributions to these series was examined. The results showed that the POT series is more suitable for the analysis of the aforementioned variables. Also, among the different distributions, the generalized extreme value (GEV) distribution has a better fitness. Next, a trend analysis and bivariate frequency analysis in two stationary and non-stationary modes was performed using actual precipitation data, maximum and minimum temperatures, and maximum sea level in the coastal area of Manhattan (USA), and the impact of global warming on coastal flooding was investigated. Also, stationary copula functions were used to examine the dependence structure, estimate the bivariate density distribution of sea level with precipitation and temperature. After determining the appropriate copula function, flood values were determined with a joint return period of 2 to 1000 years and compared with observational events in the past. The results of trend analysis indicated a significant increasing trend in sea level data, and a decreasing trend in the interval time between occurrences of joint events. baxsed on the results of the univariate stationary frequency analysis, the 1000-year sea level equivalent to Hurricane Sandy was estimated to be 5.52 m, while the multivariate frequency analysis showed that the maximum sea water level rise, combining temperature and precipitation, was 1.5 m. Also, in this thesis, a method for non-stationary multivariate frequency analysis of the maximum sea level rise was developed using hydroclimatic variables as auxiliary variables of the univariate generalized extreme value (GEV) distribution, as a marginal distribution function and copula functions. The results showed that for the maximum sea level rise, the location parameter of the marginal distribution is directly related to the maximum temperature auxiliary variable. Also, for precipitation, the scale parameter is related to the minimum temperature auxiliary variable, and the shape parameter is time-dependent. For this purpose, a non-stationary marginal distribution was created with generalized additive models (GAMs) for location, scale, and shape parameters (GAMLSS). The univariate return periods of Hurricanes Sandy and Irene in the non-stationary marginal distribution of GEV were estimated to be 85 and 12 years, respectively, while for the stationary marginal distribution of GEV, these values were estimated to be 1200 and 25 years, respectively. In order to analyze the joint frequency of precipitation and sea water level, the non-stationary marginal distribution fitted to these variables was placed in the time-varying copula function. The results showed that in the bivariate frequency analysis of maximum sea level rise and precipitation, the normal copula function has more flexibility than other copula functions. Using the time-varying copula function, the bivariate return periods of Hurricanes Sandy and Irene were obtained to be 109 years and 136 years, respectively. Next, in order to analyze the regional frequency, by combining the data of the 8 selected stations, the generalized gamma distribution was fitted to the time series of the variables under study. After fitting the generalized gamma distribution function with time-dependent parameters to the variables of precipitation, sea water level, duration, and the interval time between the coincidences, a regional multivariate distribution was created using regular tree copula functions (C-vine). Then, using the developed non-stationary multivariate distribution, the multivariate return period of regional extreme events was calculated. These results indicate that the analysis of the joint frequency of precipitation and maximum sea level rise in the non-stationary state can increase the accuracy of flood estimation in coastal areas, which will provide useful information to designers and planners to reduce flood losses in these areas. Although the proposed method can also be applied to other hydroclimatic variables, the findings of this study indicate the necessity of considering non-stationarity in the frequency analysis of hydrological extreme events.