Hydroclimatic Trends and Streamflow Response to Recent Climate Change: An Application of Discrete Wavelet Transform and Hydrological Modeling in the Passaic River Basin, New Jersey, USA

Abstract

The exigency of the current climate crisis demands a more comprehensive approach to addressing location-specific climate impacts. In the Passaic River Basin (PRB), two bodies of research—hydroclimatic trend detection and hydrological modeling—have been conducted with the aim of revealing the basin’s hydroclimate patterns as well as the hydrologic response to recent climate change. In a rather novel application of the wavelet transform tool, we sidelined the frequently used Mann–Kendal (MK) trend test, to identify the hidden monotonic trends in the inherently noisy hydroclimatic data. By this approach, the use of MK trend test directly on the raw data, whose results are almost always ambiguous and statistically insignificant in respect of precipitation data, for instance, no longer poses a challenge to the reliability of trend results. Our results showed that, whereas trends in temperature and precipitation are increasing in the PRB, streamflow trends are decreasing. Based on results from the hydrological modeling, streamflow is more sensitive to actual evapotranspiration (ET) than it is to precipitation. In periods spanning decades with sufficient water availability, energy governs actual evapotranspiration rates, rendering streamflow more sensitive to increases in precipitation. Conversely, during meteorologically stressed decades, water availability dictates actual evapotranspiration, consequently amplifying streamflow sensitivity to fluctuations in actual evapotranspiration. We found that the choice of baseline condition constitutes an important source of uncertainty in the sensitivities of streamflow to precipitation and evapotranspiration changes and should routinely be considered in any climate impact assessment.

1. Introduction

Global climate change is expected to accelerate the global hydrologic cycle, which will drive more intense floods and droughts leading to changes in streamflow and water resource availability. An alteration of the discharge regime of rivers [1,2] is usually the ultimate consequence. In the past decades, empirical evidence of warming-driven intensification of the hydrologic cycle has led to an increasing interest in the linkage of climatic variability or change to hydrological processes across space and time [3]. More often than not, the literature is either rich in the detection and analysis of hydroclimatic trends (e.g., [4,5,6]) or hydrological modeling studies (e.g., [7,8,9]), without considering both. For instance, ref. [10] underscored the importance of understanding climate variability and trends for the management and planning of water resources. Ref. [11] also investigated monthly and annual trends in temperature and precipitation and evaluated the significance of their variability for crop yields using multiple regression analysis. In another study, ref. [12] used the Soil and Water Assessment Tool (SWAT) model combined with land use and global climate models to examine how climate and land use/land cover work together to impact blue-green water in an arid basin. At a time when the global warming problem has evolved into a crisis [13,14], it is important that hydrological impact assessments be carried out from a holistic standpoint. Although the analysis and detection of trends can provide useful insights in terms of a general estimate of the direction and changes in magnitudes of hydrometeorological series, they lack the ability to predict unprecedented future conditions. Process-based models, although only capable of representing processes to the scope that they are quantitatively understood, can provide a robust framework for assessing hydrological response to climate change [15]. More so, because the direction and extent of changes in river flows are dependent on the relative balance between precipitation and the processes that govern evapotranspiration [15], the causes of discharge changes—which oftentimes seem controversial [16,17]—can effectively be examined from a hydrological modeling standpoint.

In studies that emphasize the detection and analysis of hydroclimatic trends (e.g., [16,17,18,19]), one statistical tool that has commonly been used is the nonparametric Mann–Kendall (MK) test and its modified forms. The MK trend test identifies changes in hydroclimatic series by simply fitting a monotonic (e.g., linear) trend at a certain time-period where a significant level is assigned by a statistical test. While the robustness of this test is not in doubt, its application to hydroclimatic time-series can be particularly challenging due to the nonmonotonic and nonuniform character of hydroclimatic variables. In many of these studies, trends in hydroclimatic variables, especially, precipitation, is found to be either not statistically significant or significant at a lower confidence level. For instance, in studying trends in precipitation, temperature, and streamflow at 13, 12, and 9 gauging stations, respectively, ref. [20] found that mean annual and seasonal precipitation trends were not statistically significant in all the gauging stations. In the entire Swat River Basin of Pakistan, ref. [21] observed that annual precipitation time-series did not show any statistically significant trends in all the subbasins examined. Ref. [6] also found that no statistically significant monotonic trend was detected for annual rainfall, although nonsignificant downward trends were dominant. Numerous studies with similar findings abound in the literature, and are usually in the application of Mann–Kendal tests to hydroclimatic time-series (e.g., [22,23,24,25,26]). Ref. [27] noted that, because the stochastic structure of time-series data has the tendency to assume trend-like features, analyzing trends in nonstationary time-series can cause a purely stochastic behavior to appear deterministic, leading to a likely erroneous interpretation of results. More so, because climatic phenomena and events (e.g., precipitation, hurricanes) are products of various complex atmospheric processes [28], noise is inevitably present, and this can affect the variability and trend in the data series. In hydroclimatic times series where nonmonotonicity is more the rule rather than the exception [29], identifying the hidden monotonic trend and assessing their statistical significance subsequently provide more reliable results than those derived from the direct trend analysis of the raw data [30].

Because the structure of hydroclimatic data is often hidden behind the noise, a precise mathematical operation that looks at the data through the noise and quantifies the structure present in the signal is needed. One such tool is the wavelet transform (WT). WT is a relatively recent development in the field of signal processing [31,32], and has, in recent times, emerged as an effective tool to analyze trends in hydroclimatic series especially in the atmospheric and hydrological science space [29,33,34,35,36,37,38]. It can be thought of as a ‘mathematical microscope’ with the ability to zoom in and out of the signal (or time series) to pull out the patterns. In its application, a signal or time-series data are decomposed into their low-frequency components and high-frequency components. The different decomposition levels, representing different periodic time scales, are subsequently analyzed for trends. The last decomposition level, which contains the lowest frequency component, usually represents the trend component of the time series. Thus, among the methods presently used in analyzing time-series data, the wavelet approach has the superior ability to handle the nonstationary characteristics of hydroclimatic time-series on multiple temporal resolutions [39], making it well suited for identifying trends over a long period of time.

In light of understanding and quantifying the hydrological impacts of climate variability/change, different approaches (i.e., conceptual methods, analytical, experimental, and hydrological modeling) have been used. Among these, process-based hydrological models provide a means to examine the physical mechanisms and processes that drive hydrological changes and variations. Their primary purpose is to partition precipitation into evapotranspiration and streamflow. They must, however, be thoroughly evaluated against field observations that sufficiently represent the region and timeframe of interest [40,41]. By applying a rigorously calibrated and validated physically based hydrological model, MIKE SHE, to the Rockaway catchment, a subbasin of the Passaic River Basin (PRB) in New Jersey, USA, we explored the mechanisms underpinning streamflow changes through the examination of MIKE SHE simulated water balance terms under various climate scenarios.

Thus, the novelty of this study lies in the application of advanced trend analysis tool with a physically based hydrological model that simulates both surface and subsurface flows in the land phase of the hydrological cycle. This combination will provide important clues on the key underlying variables behind the trend as well as insights into how hydroclimatic patterns may change into the future.

In the PRB and its surrounding areas, the lingering effects of a troubled history of improper environmental practices from the industrial boom continue to be experienced. According to [42], the Passaic River played a central role in the early development of New Jersey. In the late 18th century, the river served as navigable routes connected by a system of canals to the Delaware River. It was also an early source of hydroelectric power at the Great Falls in Paterson, making the region a focal point for industrial mills. Consequently, the lower Passaic experienced significant environmental contamination resulting from decades of industrial activity in the vicinity. By 1970, issues pertaining to flooding had already become apparent due to the presence of dams, and they continue to afflict the inhabitants of the basin at present. Furthermore, the intricate network of river systems amidst the heterogeneous biophysical landscape within the basin presents a complex array of conflicting interests and water-related challenges. In a region where increases in temperature [43,44] accompanied by evapotranspiration and snowpack depletion [15,45] are both observed and projected, the consequent impact on streamflow could be extensive. In light of these circumstances, our study undertakes a comprehensive hydrological impacts assessment in pursuit of establishing important foundations for the predictive understanding of the impacts of climate change on water resources in the PRB and its environs.

To this end, this research endeavor embraces a dual-fold objective: (1) Identifying shifts in hydrometeorological trends within the PRB spanning the period from 1979 to 2021; and (2) Investigating the responsiveness and sensitivities of hydrological systems to climatic phenomena in the Rockaway River basin, a subcatchment of the PRB. The analysis of hydroclimatic trends employs the discrete wavelet transform and Mann–Kendal test methodologies, while the evaluation of climate change implications in the PRB employs a hydrological modeling approach coupled with sensitivity analysis.

The remaining part of the paper is organized as follows. Section 2 describes the study area and data source. Section 3 details the methodology used, which includes the discrete wavelet transform, hydrological model evaluation, and hydrologic impacts assessment. Results of the trend analysis and hydrological model performance, and model assessment are presented under results and discussion in Section 4 followed by the conclusion in Section 5.

2. Study Area and Data Source

2.1. Study Area

The nontidal portion of the PRB is elliptical in shape, draining approximately 2135 square kilometers of Northern New Jersey (NJ) and Southern New York State (NY). It is bounded by longitude 74°1′1″ and 74°39′16″ W and latitude 40°35′23″ and 41°23′37″ N, intersecting six (6) counties in NJ, two (2) in NY (Figure 1). The entire basin stretches across three (3) Watershed Management Areas (WMA-03-04-06) with seven major tributaries: Whippany River, Rockaway River, Pompton River, Pequannock River, Wanaque River, Ramapo River, and Saddle River. Physiographically, the basin can be divided into three main regions: the series of parallel ridges that trend northeast/southwest forming the Highlands; the Central Basin, comprising large areas of swamps and meadows; and the roughly flat Lower Valley. Winding through seven counties and 45 municipalities, the Passaic River originates from near the Borough of Mendham (Morris County), and finally empties into the Newark Bay. The nontidal part of the river is regulated by 10 major reservoirs (Canister, Greenwood lake, Clinton, Oak Ridge, Charlotte-burg, Echo lake, Split Rock, Monksville, Wanaque, and Boonton) to provide flood control and water supply, among other things, to surrounding municipalities. Aggregated reservoir storage in the basin is about 68,533 million gallons (MG). Of the reservoirs, Wanaque reservoir is the largest, with storage capacity of 29,630 MG (43%) of the basin total (New Jersey Water Science Center, 2018; Survey, 1964, 1970; Wells, 1960). A summary of the hydrometeorological conditions in the basin is given in

Table 1. Basic hydroclimatic information in the Passaic River Basin.

	Drainage Area (sqkm)	Area (% of PRB)	Annual Flow (m3s)	Temperature (°C)	Precipitation (mm)
			Mean (Min–Max)
PRB	2135	-	402,088 (30,968–958,992)	10.59	1281
RA	300.4	14.07	2513 (611–4037)	9.52	1296
RW	46.4	2.17	337 (122–721)	9.74	1298
UP	356.3	16.67	1916 (344–2977)	11.11	1269

PRB (1983–2021 WY); RA (1971–2010 WY); RW (1986–2021 WY); UP (1971–2010 WY).

. Mean annual flow at the outlet of the basin is estimated at 402,088 m³/s for the 1983–2021 period. On average, the Rockaway river (RA) contributes about 0.62 percent of flow to the Passaic river. The study basin lies within the modified continental climate zone, characterized by hot summer and cold winter (Paulson, 1991). Moving from north to south in New Jersey, the modified climate zone comprises five (5) main divisions: North, Central, Southwest, Pine Barrens, and Coastal zones with PRB located in the North and Central climate zones. For the period 1981–2010, mean annual precipitation of 1281 mm (50.4 in) occurred over the PRB with higher values (1298 mm or 51.1 in) in the Ringwood catchment and the lower values (1269 mm or 49.96 in) in the Upper Passaic (UP) catchment. Mean temperature for same period in the study basin is calculated as 10.59 °C. Colder temperatures are observed over the RW and RA catchments whereas hotter temperatures occur in the UP area. Throughout the PRB, mean annual actual evapotranspiration is estimated to be approximately 793 mm or 31.2 in [46].

2.2. Hydrometeorological Data

For this study, the widely used gridded observations from Parameter-elevation Regression on Independent Slopes Model (PRISM, Oregon State University, http://prism.oregonstate.edu (accessed on 19 September 2022)) provided meteorological data. Flow data for the Rockaway and Upper Passaic subcatchments were obtained from records of reconstructed streamflow by [47] whereas data for the Ringwood catchment were sourced from the United States Geological Survey (USGS) water data website (USGS, https://waterdata.usgs.gov/nwis/ (accessed on 25 November 2022)).

The data used in this study spanned the periods 1979–2021 and 1981–2021 water years (WY) for the trend analysis and hydrological modelling, respectively. Missing flow data, when present, were handled based on streamflow outputs from a duly calibrated and validated hydrological model of the subcatchment [Ringwood, correlation coefficient: 0.85 and Nash–Sutcliffe: 0.71]. In all, seven (7) hydrometeorological variables (i.e., flow, m³/year, precipitation [precip, mm/year], minimum temperature [Tmin, °C], mean temperature [Tmean, °C], maximum temperature [Tmax, °C], number of days with precipitation greater that 10 mm [R10, day], and consecutive dry days [CDD, day]) for the three (3) studied subcatchments were processed and aggregated into annual time scales for the trend analysis.

2.3. Land-Use, Soil and Elevation Data

In hydrological processes, the combined effect of land cover, soil, elevation, and other catchment characteristics are reflected in the flow dynamics of river systems in a basin. The 2011 Land cover data were available from the National Land Cover Dataset (NLCD) [https://www.mrlc.gov/viewer/ (accessed on 13 July 2022)], and simplified into six (6) dominant land cover/vegetation classes (i.e., developed, forest, agricultural, wetlands, bare land, and water). Soil information was accessed using the United States Department of Agriculture (USDA) soil data viewer software [version 6.2]. The topography of the PRB was defined by a digital elevation model (DEM) extracted from the USGS database at 10 m spatial resolution (https://apps.nationalmap.gov/downloader/ (accessed on 15 May 2020)).

3. Methods

In line with the objectives of the study, two major tasks were carried out: 1) the analysis of hydroclimatic trends via the discrete wavelet transform (DWT) approach; and 2) the development of a hydrological model to assess the impacts of recent climate changes on water balance terms (i.e., precipitation, evapotranspiration, and streamflow). The hydroclimatic trend detection was conducted for seven (7) indicator variables in three (3) subcatchments because of their physiographically distinct locations. The hydrological modeling study was, however, conducted for only the Rockaway sub-basin given its relatively large size, largely representative characteristics to the PRB, and available groundwater data. Summary of the steps involved in our analyses are outlined below and described in more detail in the forthcoming sub-sections:

• Seven (7) different hydroclimatic indicator variables used in the trend analysis were derived from temperature, precipitation, and streamflow data series obtained for each subcatchment. They were mean annual Tmin, Tmean, Tmax, Precip, Flow, R10, and CDD spanning the period 1979–2021.
• Each time series was decomposed via the DWT, having selected the Daubechies (db) wavelet, deemed an appropriate mother wavelet in our study context, to split the series into their high frequency detailed (D) and low frequency approximate (A) components.
• The MK Z-values of the original signal and the approximation of each Daubechies (db) wavelet form starting from db4–db10 (e.g., [48]) were computed to determine the wavelet form that gives MK Z-value closer to that of the original signal. This was the optimal trend from the approximation components of each analyzed time series.
• Having selected the optimal monotonic trend, an MK test was subsequently applied to determine the statistical significance of the DWT-based trend.
• A hydrological model for the Rockaway subbasin was developed, calibrated, and validated, and the performance of the model against the observed streamflow and groundwater data was evaluated using standard statistical criterion. The water balance module was run to obtain outputs of water-balance components for the impacts assessment.
• Change point analysis was carried out to divide data into the naturalized or baseline periods, where minimum effects of human activity on streamflow is expected and impacted periods. Subsequently, a climate elasticity exercise was undertaken to explore sensitivities of climate variables to streamflow and corresponding contributions in the Rockaway sub-basin.

3.1. Discrete Wavelet Transform

Wavelet transform (WT) is a mathematical tool that uses wave functions known as wavelets—akin to the sine and cosine functions in Fourier transforms (FT), to convert signals or time-series data into different frequency components. WT rides on the fundamental concept of Fourier transform, which operates on the idea that any function can be decomposed into a sum of pure waves with different frequencies. Therefore, the frequency domain represents the relative contributions of each frequency that comprises the function. The major limitation with FT is that knowledge about frequency is accessed at the expense of the temporal dynamics (i.e., there is no clue as to when certain frequencies begin or end). As a result, wavelets come into play to resolve this inherent trade-off of information between frequency and time in the FT [49]. Through the application of wavelet transform, an optimal frequency–time balance is attained. The key feature about wavelets is that the wave-like oscillations are short-lived and localized in time. It is worth noting that a wavelet is not just a function, but a whole family of functions which all satisfy certain requirements. The popular family of functions include Daubechies, Coiflet, Symlet, Haar, Morlet, Gaussian, Shannon, Meyer, and Mexican Hat; each one of these is tuned for specific applications. In general, to be considered a proper wavelet, a function must satisfy two main constraints: (1) the admissibility condition of having a zero mean, and (2) the finite energy condition of having a limited duration, from which a function attains its localized nature in time. In short, wavelet analysis is a completely flexible windowing technique that allows a function to change over time based on the shape and compactness of the time-series signal [50]. By this very nature, different modes of variability that varies in time can be extracted in the WT process, allowing the time-frequency characteristics of any kind of signal to be analyzed (Wei et al., 2012). Recent years have seen a wide range of studies using WT, especially to analyze hydrometeorological time-series (e.g., [29,30,51,52,53,54,55]). While the vast majority of these studies focused on trends, others emphasized the dominant periodic time scale responsible for the trends.

In the WT process, as a mother wavelet moves across the signal, several coefficients are generated according to the similarity between the signal and the mother wavelet at any specific scale. Generally, WT is divided into two main types, the continuous wavelet transform (CWT) and the discrete wavelet transform (DWT). The continuous type can generate quite numerous and often redundant coefficients at every resolution level, making its application and interpretation more complex and uncertain. The DWT is, however, considered a more effective approach, having the ability to overcome the data redundancy issue by simplifying the transformation process based on the dyadic (power of 2) scale [56]. Given a suitable wavelet family and decomposition level, the DWT decomposes a series into several sub-series during the transformation process (Whitcher et al., 2002). Following Equation (1) below (reader is referred to [56] for details), the coefficients of DWT can be calculated:

$$W_{\phi }\left(a,b\right)=\frac{1}{{\left(2\right)}^{a/2}}{{\displaystyle \sum }}_{t=0}^{N-1}X\left(t\right)\phi \left(\frac{t}{2^a}-b\right)$$

(1)

where $2^a$ denotes the dyadic scale of the DWT. Note that the resulting detail and approximation coefficients from the decomposition are merely intermediate coefficients, and has to be reconstructed, first to their approximation and detail components and then to the original signal. This readjustment to the original one-dimensional signal ensures that each component has the same length as the original signal, thereby enabling proper investigations of their contribution to the signal [57,58,59]. In a simplified form, the reconstruction of the detail and approximation components can be computed as:

$$S\left(t\right)=A_n\left(t\right)+{{\displaystyle \sum }}_{l=1}^n{D}_l\left(t\right)$$

(2)

where $S\left(t\right)$ is the original signal, and $A_n\left(t\right)$ is the approximation component at level $n$, and $D_l\left(t\right)$ is the details component at different levels (where $l$ = 1, 2, 3, …, $n$ denotes index for the levels). In MATLAB, computation of a perfect signal reconstruction is achieved using the Inverse Discrete Wavelet Transform (IDWT).

Although nearly all hydroclimatic processes are continuous in nature, their available time-series outputs are delivered in discrete formats (Wilks, 2011), making its use with DWT more appropriate than that of CWT. In the application of DWT, the original time-series signal is passed through low-pass and high-pass filters and emerge as Approximation (A) and Detail (D) components, respectively. While component D represents the small scale, high-frequency series, component A comprises the high scale, low-frequency series [34]. The decomposition process can continue iteratively, where component A from the first decomposition is further divided into new A and D components [33,49,54,60]. In this study, the Daubechies mother wavelet was chosen because of its characteristic orthogonality, and compact support, which are very important properties for localizing events in signal analysis, and deemed appropriate for hydrometeorological time-series [33,61].

Time-Series Decomposition via DWT

In our wavelet analysis, the one-dimensional flow signal (data series) and each of the temperature and precipitation indicator variables served as inputs to the multi-level 1-D wavelet decomposition function in the MATLAB Wavelet Toolbox (MATLAB version R2021a). With the db wavelet family as the mother wavelet, the operation produced a wavelet transform of each input time-series signal at all dyadic scales. Three main parameters were taken into account during the DWT process [33]: (1) the appropriate type of db wavelet; (2) a suitable signal border extension method; and (3) the most appropriate number of decomposition levels. First, several forms of db wavelet (e.g., db1–db10) exist, and the appropriate type must be selected for the decomposition process [62]. As suggested by [33], a useful method in selecting the appropriate db wavelet type is to calculate the relative error (R_E) between the MK Z-values of the original signal and that of the approximation (A) of the last decomposition level. R_E is computed as follows:

$$R_E=\frac{\left|Z_a-Z_o\right|}{\left|Z_0\right|}$$

(3)

where $Z_a$ and $Z_o$ are the MK Z-values of the approximation of the last decomposition level and the original dataset, respectively. For each indicator variable and study catchment, the appropriate db wavelet was selected to minimize $R_E$. Because trends are supposed to be gradual and slowly changing process, smoother db wavelets (i.e., db4–db10), considered as better in detecting time varying behavior over the long term [33,63], was used in the selection of the appropriate db wavelet in the study (

Table 2. Daubechies (db) wavelet type, minimum relative error (RE), Mann–Kendal test, and Sen’s slope (SS) for each metric in each subcatchment (* denotes significant at p = 0.05).

Subcatchment	Parameter	Metrics
		Precip	Flow	R10	CDD	Tmin	Tmean	Tmax
Ringwood	Wavelets	db7	db6	db8	db5	db4	db4	db4
	RE	3.75	0.03	4.86	5.72	3.93	4.48	13.46
	MKSL	871 *	−457 *	−457 *	−877 *	903 *	903 *	903 *
	SS	0.723	−0.165	−0.024	−0.051	0.047	0.033	0.018
Rockaway	Wavelets	db5	db6	db8	db4	db4	db4	db10
	RE	7.29	4.18	2.91	3.02	4.43	4.29	10.33
	MKSL	635 *	−745 *	577 *	293 *	903 *	903 *	831 *
	SS	0.129	−2.406	0.083	0.0042	0.059	0.035	0.0034
Upper Passaic	Wavelets	db4	db4	db7	db4	db4	db4	db4
	RE	7.49	6.32	67.10	7.79	3.21	3.79	6.86
	MKSL	903 *	903 *	433 *	213 *	903 *	903 *	903 *
	SS	2.401	7.712	0.0134	0.0062	0.0375	0.0253	0.013

). Second, border extension is an important consideration due to the issue of border distortions in the DWT process, arising because of the finite length of the signal. Thus, the decomposition process cannot occur outside the two limits (i.e., the start and end points) of a signal as there is no available information beyond the ends [64]. As suggested by [64], three different border extension methods are employed to address the issue: zero-padding, symmetrization, and periodic padding. In our analysis, symmetrization—which is the default mode in MATLAB, was used. It assumes that signals beyond the original support can be retrieved by symmetric boundary replication [54]. Finally, the relevant number of decomposition level must be determined in order to avoid unnecessary levels of data decomposition especially, for larger datasets (see [33,54]). This will, however, depend on the length of data points as well as the type of mother wavelet used. According to [65], the maximum number of decomposition level, $L$, can be calculated from Equation (4) below.

$$L=\frac{Log\left(\frac{n}{2v-1}\right)}{Log\left(2\right)}$$

(4)

where $n$ is the number or length of data points in the time series and $v$ is the number of vanishing moment of a db wavelet [34]. In MATLAB, the number of vanishing moments ($v$) is equal to the db wavelet type number (i.e., 1–10). Note that the number of data points ($n$) in a time series is not exactly in a dyadic format (as in the case of this study, 43 data points). Thus, the DWT computation in MATLAB is carried out using the nearest upper dyadic arrangement. Therefore, the maximum decomposition level based on our data points was calculated to be 6 in the study. Additionally, because data decomposition via DWT assumes a dyadic format, each of the decomposed component represents a different period of integer powers of two from the lowest scale. Therefore, D1, D2, and D3, respectively, represents 2-, 4-, and 8-unit periodic component, in that order, according to the time scale (e.g., seasonal, monthly, annual) used in the analysis. For example, D2 will represents a 4-year or 4-month intervals in an annual or monthly data series, respectively, but 12-month intervals for a seasonal data series since its time step is 3 months.

3.2. Trend and Change-Point Detection Tests

In this study, the Mann–Kendall (MK) test [26,35,66,67,68] was used to determine the statistical significance of the DWT-based trend [48]. It is probably the most widely used nonparametric statistical test for monotonic trend evaluation [69]; and well noted for its simplicity, robustness, and resilience to missing values in a data series [63]. One key issue that may arise when using the MK test is the presence of serial correlation or autocorrelation, very common in precipitation and streamflow data [56]. It occurs when a variable and a lagged version of itself is observed to be correlated between two successive time intervals. If the lag-1 autocorrelation in a time series is found to be significant, the modified MK test must be used [70]. Although autocorrelation issues are not common in annual data series, we applied the modified version of the MK test in this study where significant lag-1 autocorrelation was detected in our data series.

Furthermore, a change-point analysis was performed to identify the most likely year(s) in our streamflow data where significant changes could occur [71,72]. This was key in our hydrological impacts analysis where we needed to explore naturalized periods when stream flow experienced little or no disturbance as well as impacted periods. Various change-point methods exist, including the sequential Mann–Kendal test [73], Pettit’s test [74], the cumulative sum (CUSUM) test [75], and the Worsley Likelihood Ratio Test [76]. Using the R packages changepoint and ecp [77], the distribution free CUSUM test [78] and the complementary Permutation test [79] were used because they revealed similar break-point years in our streamflow time-series. These methods detect significant changes in the mean or distribution of a time series when the exact times of the changes are not known. For a detailed description of these methods, we refer the reader to the relevant literature cited.

3.3. Hydrological Model Development for the Rockaway River Basin

In the present study, the numerical code used for our hydrological assessment studies is MIKE SHE [80]. MIKE SHE is an integrated, fully distributed, physically based hydrological modeling system [81,82], that simulates all the major hydrologic process in the land phase of the hydrological cycle including evapotranspiration, overland flow, unsaturated flow, saturated flow, and streamflow (Figure 2). It uses the hydrodynamic model MIKE 11/MIKE Hydro to simulate channel flow and lakes (using flood code) in one dimension. For a detailed description on the development and modelling structure of the MIKE SHE hydrologic model the reader is referred to the MIKE Zero user manual by [81].

Model Calibration and Validation

A distributed hydrologic model such as MIKE SHE typically requires large number of model parameters to be assigned. Although these parameters have a clear physical meaning and can be defined explicitly from field measurements, ref. [84] suggested that the number of parameters subject to calibration should be as small as possible. For this study, initial values as well as ranges of primary parameters from field data, the published literature, and prior modelling experience guided the calibration process. The manual “trial and error” procedure was first applied, which involved perturbing one parameter while keeping all other parameters unchanged. This was carried out repeatedly within a reasonable range of values for a series of model runs until a favorable agreement between measured and simulated flow and groundwater level was achieved. Following the manual approach, an automatic calibration was conducted. Finally, validation was carried out to ensure that model parameters derived from calibration were generally valid.

Prior to the model calibration, change point analysis was performed on the streamflow data for the time span of 1981–2022 to find likely break-point year(s). Accordingly, the data were divided into baseline periods (1982–1991, 1992–2001, and 1982–2001) and impacted periods (2002–2011 and 2012–2021), as mentioned earlier. The calibration and validation of the model was carried out within the baseline period for 1982–1986 and 1986–1991, respectively. Typically, the simulation period includes the first few months of warm-up period to stabilize the model; as well as the calibration and validation periods. The adequacy of the model was evaluated based on four standard statistical criteria used in MIKE SHE: mean error (ME), root mean square error (RMSE), correlation coefficient (R), and the widely used Nash–Sutcliffe coefficient (NSE). These indicators detect system errors and the goodness of fit between simulated and observed monitoring observations in the form:

$$M{E}_i=\frac{{{\displaystyle \sum }}_t\left(Ob{s}_{i,t}-Cal{c}_{i,t}\right)}{n}$$

(5)

$$RMSE=\frac{\sqrt{{{\displaystyle \sum }}_t{\left(Ob{s}_{i,t}-Cal{c}_{i,t}\right)}^2}}{n}$$

(6)

$$R=\sqrt{\frac{{{\displaystyle \sum }}_t{\left(Cal{c}_{i,t}-{\overline{Obs}}_i\right)}^2}{{{\displaystyle \sum }}_t{\left(Ob{s}_{i,t}-{\overline{Obs}}_i\right)}^2}}$$

(7)

$$NSE=1-\frac{{{\displaystyle \sum }}_t{\left(Ob{s}_{i,t}-Cal{c}_{i,t}\right)}^2}{{{\displaystyle \sum }}_t{\left(Ob{s}_{i,t}-{\overline{Obs}}_i\right)}^2}$$

(8)

where t is the simulation time in day; n is the total simulation days; i is the calibration point i; Obs_i,t is the observed daily discharge at location i at day t; $\overline{Ob{s}_i}$ is the mean of the observed discharge at location i for the simulation period, and Calc_i,t is the simulated discharge at location i at day t.

3.4. Hydrological Impacts Assessment

After successfully calibrating and validating the hydrologic model for the period (1982–1991) considered to be within the naturalized undisturbed periods, the model was run with climatic inputs to simulate discharge for both the naturalized periods and the impacted periods identified by the change point analysis. In all, discharge for five (5) different periods were simulated, and a water balance output obtained for precipitation (P), streamflow (Q), and actual evapotranspiration (ET). Further, we assessed the hydrological impacts by computing changes between the baseline periods and impacted periods for the water balance components. Finally, the concept of elasticity as proposed by [85] was employed to evaluate the sensitivities of streamflow to changes in climate. According to this concept, climate elasticity of streamflow is the proportional change in streamflow divided by the proportional change in a climate variable. For instance, the precipitation elasticity of streamflow is defined as:

$${\epsilon }_p=\frac{dQ/Q}{dP/P}=\frac{dQ}{dP}·\frac{P}{Q}$$

(9)

Likewise, the actual evapotranspiration elasticity of streamflow is:

$${\epsilon }_{ET}=\frac{dQ/Q}{dET/ET}=\frac{dQ}{dET}·\frac{ET}{Q}$$

(10)

In applying Equations (9) and (10) to the water balance outputs obtained for the Rockaway catchment model, the relative contributions of precipitation and actual ET changes to streamflow changes can be quantified. Over the long term, the water balance model can be expressed as [86]:

$$\overline{Q}=\overline{P}-\overline{ET}$$

(11)

where $\overline{Q}$, $\overline{P}$, and $\overline{ET}$ denote long term mean values. In Equation (11), there is an implicit assumption that groundwater flow into and out of the Rockaway subcatchment cancels out and storage change over the long term is negligible.

For a largely undisturbed catchment, the changes in streamflow between two periods ($dQ$) based on Equation (11) can be estimated as:

$$dQ=d{Q}_P+d{Q}_{ET}$$

(12)

with $d{Q}_P$ and $d{Q}_{ET}$ denoting the contribution to streamflow from precipitation and actual ET, respectively. Combining Equations (9), (10) and (12), $dQ$ can be rewritten as:

$$dQ=d{Q}_P+d{Q}_{ET}=\left({\epsilon }_{P·dP/P}+ {\epsilon }_{ET·dET/ET}\right)Q$$

(13)

where $dP$ and $dET$ are changes in precipitation and actual evapotranspiration between two periods. ${\epsilon }_p$ and ${\epsilon }_{ET}$ are precipitation elasticity and actual evapotranspiration elasticity of streamflow, respectively. According to [87], a 1% change in P or ET triggers an ${\epsilon }_p$ or ${\epsilon }_{ET}$ percent change in Q.

Note that flow data used in this analysis form part of the reconstructed streamflow records by [47] for selected watersheds in the PRB, and therefore the Rockaway subcatchment is assumed to be largely undisturbed for the purpose of this study. Being mindful of the fact that climate elasticity to streamflow varies depending on the location and reference period [88], we explored the sensitivities of flow relative to different baseline periods. Thus, Equation (13) was set up for the two impacted periods (2002–2011 and 2012–2021) relative to three baseline periods (i.e., 1982–1991, 1992–2001, and 1982–2001). The values of ${\epsilon }_p$ and ${\epsilon }_{ET}$ were then computed simultaneously from two equations to obtain the contributions of P and ET changes to streamflow change in the Rockaway subbasin.

4. Results and Discussion

4.1. Decomposition of Time-Series Data via DWT

According to Equation (3), the types of db mother wavelet that produced the optimal parameters for the decomposition process for each dataset are presented in Table 2. Figure 3 illustrates an example of the decomposition results for the flow data in the Rockaway subcatchment. The original time series or signal (S) can be reconstructed by summation of all the detailed components (D1–D6) and the approximation component of the last decomposition level (A6). It can be seen that at higher decomposition levels, the frequency of the D components decreases. The last decomposition level of the A component (A6) shows the trend of streamflow in the Rockaway catchment. On a dyadic scale, D1 depicts the time series of a 2-year mode, D2 shows a 4-year mode, D3 is in an 8-year mode, D4, a 16-year mode, D5 corresponds to a 32-year mode, and 64-year mode for D6. These modes are the time scales at which those cycles are revealed, implying that for a dataset spanning a period of 42 years, the trend as revealed by the DWT for stream flow in the Rockaway catchment could only emerge over a 64-year cycle. It was thus impossible to see this trend just by applying the MK trend test on the raw dataset. The same process was replicated in all three subcatchments for each hydroclimatic indicator variable, as depicted by Figure 4, Figure 5 and Figure 6. At a confidence level of 95% (i.e., p-value = 0.05), MK statistics were subsequently applied to the decomposed times series.

4.1.1. DWT Trend Analysis of Hydroclimatic Indicators

Results from the trend analysis using the discrete wavelet transform for Precip, Flow, R10, CDD, Tmin, Tmean, and Tmax from the Ringwood, Rockaway, and Upper Passaic catchments are shown in Figure 4, Figure 5 and Figure 6. Mann–Kendall statistics (i.e., significant level (SL) and Sen’s slope (SS)) applied on the DWT trend results are also summarized in Table 2. The positive and negative MK values indicate significantly increasing and decreasing trends, respectively, and the magnitude of the trends is described by the SS values. All the analyzed hydroclimatic signals were significant at p = 0.05, identified by the asterisk.

4.1.2. Hydroclimatic Trends in the Ringwood, Rockaway, and Upper Passaic Subcatchments

DWT trend results for all hydroclimatic variables in the Ringwood subcatchment are shown in Figure 4. Precipitation shows a significant increasing signal at a rate of 0.723 mm/year. However, this increase does not reflect in the streamflow trend in the Ringwood subcatchment. Flow is rather showing a significant downward trend beginning from 1996 through to 2021 at a rate of 0.165 m³yr⁻¹. It does appear that the downward trend observed for streamflow largely tracks with heavy precipitation (R10) rather than mean precipitation, and corroborated by the decreasing trend in consecutive dry days (CDD). The observed significantly increasing trend in minimum, mean, and maximum temperatures suggest that temperature drives the flow dynamics in the Ringwood subcatchment with minimum temperature having the highest magnitude at 0.047 °Cyr⁻¹ and maximum temperature having the least at 0.018 °Cyr⁻¹.

Long-term trends in hydroclimatic indicator variables in the Rockaway catchment are illustrated in Figure 5. Similar to the Ringwood subcatchment, precipitation and streamflow are trending in opposite directions. As precipitation trends upward, flow is trending downward at rates of 0.129 mmyr⁻¹ and 2.406 cmyr⁻¹, respectively. Quite interestingly, a significantly upward trend is observed for heavy precipitation, in line with mean precipitation, yet these increases do not reflect in the observed flow trend. Given that consecutive dry days show a significantly increasing trend in tandem with minimum, mean, and maximum temperatures, there is a likelihood that precipitation is overwhelmed by relatively high temperatures in the Rockaway subcatchment, thereby translating in the observed downward trend in streamflow.

Figure 6 shows the hydroclimatic trends in indicator variables for the Upper Passaic subcatchment. All the metrics showed significantly increasing trends over the period. The observed upward trend in mean and heavy precipitation in the same direction as flow and temperatures indicate that the hydrology of the Upper Passaic subcatchment is largely driven by precipitation rather than temperature. Precipitation and flow are increasing at a rate of 2.401 mmyr⁻¹ and 7.712 cmyr⁻¹, respectively. In the case of temperature, the rate is higher in the minimum temperature (0.038 °Cyr⁻¹), followed by the mean temperature (0.025 °Cyr⁻¹) and maximum temperature (0.013 °Cyr⁻¹) (Table 2).

4.1.3. Comparison of Hydroclimatic Trends by Catchment

Trending from north to south of the PRB, the results suggest that hydroclimatic indicator variables are spatially nonuniform in terms of magnitude and direction. Over the analyzed period (1979–2021), precipitation is observed to show increasing signals in all subcatchments. Relatively, the rate of change is observed to be rapid in the Upper Passaic subcatchment at 2.401 mmyr⁻¹ and smooth in the Rockaway subcatchment at 0.129 mmyr⁻¹. Likewise, temperatures also show significantly upward trends in all subcatchments, with mean temperature displaying the highest rate of change in the Rockaway subcatchment, followed by Ringwood, and Upper Passaic subcatchments. The observed long-term increasing trend in precipitation and temperature in the PRB is indicative of a changing climate in the basin, consistent with the dominant trends in the broader Northeast United States region (e.g., [89,90]). In terms of extremes, precipitation intensity (R10) and consecutive dry days (CDD) point towards an upward trend from north to south in the PRB, beginning from Ringwood subcatchment with a decreasing signal to increases in the Rockaway and Upper Passaic subcatchments. This observed increasing trend is also consistent with patterns in rainfall intensity in the Northeast (e.g., [91,92]), and provides further evidence to the linkage between extreme weather events and climate change.

In the case of streamflow, the results suggest that flow patterns appear to be influenced both by surface characteristics and climate in the PRB. Although trends in precipitation and temperature are observed to increase throughout the basin, the dynamics on streamflow is different, with downward trends observed in the Ringwood and Rockaway subcatchments and an upward trend seen in the Upper Passaic subcatchment. Given that Ringwood and Rockaway subcatchments lie in the mountainous heavily forested Highlands region as against the Upper Passaic in the densely populated, highly industrialized urban belt, the observed trends are not surprising. With regards to attributing causes of streamflow changes in the PRB, the sections that follow, involving the hydrological modeling study using the Rockaway subbasin as a case study, is expected to provide sufficient clues on the driving mechanism behind the flow dynamics in the study basin.

4.2. Change Point Analysis and Calibration and Validation of MIKE SHE Model

Change Point Detection

For the purpose of the hydrological impact assessment, change point detection was carried out to determine approximate years of abrupt changes in hydroclimatic time-series.

As presented in

Table 3. Estimated break point years in precipitation and streamflow in the Rockaway subcatchment (α = 0.05).

	Cumulative Sum Test	Permutation Test
Variables	Break Point	Break Point
Precipitation	1980	1979
	1990	1991
	2002	2003
	2011	2012
Streamflow	1980	1979
	1990	1991
	2002	2003
	2011	2012

, precipitation and streamflow time-series were explored in the Rockaway subcatchment using the CUSUM test as well as the permutation test, and were significant at α = 0.05. Given that river flow in the PRB is largely regulated, the similarities as revealed by the break point years for both precipitation and streamflow seem to corroborate the findings by [93] that climate change signals are apparent in both regulated and natural river systems. Following from this, the causes of streamflow were explored by examining outputs of water balance terms from the hydrological model developed for the Rockaway subcatchment. Accordingly, we performed decadal changes in water balance terms in line with the break point years given in Table 3. Over the study period (1979–2021), decades 1 and 2 spanned the periods 1982–1991 and 1992–2001, respectively, representing baseline periods 1 (BLP I) and 2 (BLP II). The overall period from 1982–2001 was also considered, denoting baseline period 3 (BLP III). For the impacted periods, 2002–2011, and 2012–2021, respectively, represents decades 3 (D III) and 4 (D IV).

4.3. Calibration and Validation of Rockaway Model

The Rockaway model was calibrated using both streamflow and groundwater level data. The full model simulation spanned the period 1981–1991 for streamflow and 2005–2016 for groundwater flow. Observed streamflow was calibrated and validated at the outlet of the Rockaway catchment along with two (2) groundwater observation wells (located in the catchment. These two wells are approximately 2.6 km apart with Berkshire Valley well having elevation of 222 m and that of Morris Maint well, 203.4 m. As shown in

Table 4. Performance criterion of calibrated and validated MIKE SHE model at the Rockaway subcatchment.

		Streamflow			Groundwater
Statistics		Calibration	Validation	Full Simulation	Berkshire Obs Well	Morris Obs Well
		1982–1986	1987–1991	1982–1991	2011–2016	2007–2012
Correlation coefficient (R)		0.85	0.87	0.85	0.83	0.28
Nash efficiency (R2)		0.72	0.71	0.72	-	-
ME		0.57	1.34	0.96	6.01	−1.86
RMSE		4.78	3.89	0.85	6.07	1.91

and Figure 7, the performance of the model was assessed using a combination of statistical indicators and graphical representation, respectively. Generally, the model can be said to have captured the evolution of the observed flow sufficiently well, with few mismatches in peak flows likely due to the gridded structure of the forcing data. Rising limbs of hydrographs and baseflow were also reasonably simulated. The resulting correlation coefficient (R), Nash–Sutcliffe (NSE), mean error (ME), and RMSE values for both calibration and validation periods are shown in Table 4. In respect of coefficient (R), Nash–Sutcliffe (NSE), the performance of the Rockaway model can be judged according to the general performance ratings as recommended by [94] (

Table 5. General range of model performance statistics [94].

Performance Indicator	Excellent	Good	Fair	Poor
Nash-coefficient (NSE)	>0.85	0.65–0.85	0.5–0.65	<0.5
Correlation coefficient (R)	>0.95	0.85–0.95	0.85–0.75	<0.75

).

In the groundwater level simulation, the transient dynamics of water level were satisfactorily simulated for the Berkshire Valley well, but poor in the Morris Maint well (Figure 8). Aside from possible errors from the DEM, the observed bias could be linked to boundary conditions at the border that favored the Berkshire Valley well more than the Morris Maint well. Although individual biases such as this are inevitable, the multiple mode of calibration (i.e., using both streamflow and groundwater data) in this study allows for simultaneous optimization of model parameters to ensure proper balance between the two solutions (i.e., simulated hydrograph and groundwater level dynamics). Thus, in general, the MIKE SHE performed reasonably well in capturing the observed streamflow and groundwater levels in the Rockaway subcatchment. On the basis of these results, we explore and quantify the possible mechanisms behind the observed streamflow changes in the sub-basin.

Relative to the three different naturalized and baseline periods (i.e., 1982–1991, 1992–2001, and 1982–2001), decadal changes in hydrometeorological variables were computed from the impacted periods (i.e., 2002–2011 and 2012–2021) (

Table 6. Mean annual changes in climatic and streamflow variables of the Rockaway subbasin for the baseline (BLP) and impacted periods (D).

Period	Tmin (°C)	Tmean (°C)	Tmax (°C)	Precip (mm)	Evapo (mm)	Flow (m3)
BLP I: 1982~1991	2.76	9.27	15.78	1306	805	2244
BLP II: 1992~2001	2.99	9.40	15.82	1208	777	1898
BLP III: 1982~2001	2.87	9.34	15.80	1257	791	2071
D III: 2002~2011	4.31	10.09	15.89	1427	826	2277
D IV: 2012~2021	5.24	10.56	15.88	1282	837	1947
D III minus BLP I	1.55	0.82	0.11	9.29%	2.58%	1.46%
D IV minus BLP I	2.49	1.29	0.11	−1.83%	3.98%	−13.22%
D III minus BLP II	1.32	0.69	0.07	18.13%	6.25%	19.94%
D IV minus BLP II	2.25	1.16	0.07	6.11%	7.71%	2.59%
D III minus BLP III	1.44	0.75	0.09	13.54%	4.38%	9.93%
D IV minus BLP III	2.37	1.23	0.09	1.99%	5.81%	−5.97%

).

4.3.1. Hydroclimatic Response to Changes Relative to 1982–1991 Baseline, BLP I

Decadal changes (i.e., 2002–2011 and 2012–2021) in Tmin, Tmean, and Tmax for the impacted periods relative to the reference showed increase in both decades, with the recent decade (D IV) being the warmest (Table 6). Compared to mean and maximum temperatures, minimum temperature is observed to have a higher increasing rate (2.49 °C/decade), and is indicative of a rapidly warming climate in the basin. Along with the increasing temperatures, actual ET is also observed to increase in both decades (2.58% and 3.98% for D III and D IV resp.). However, precipitation is observed to increase in D III but decreased in D IV. This decrease in precipitation, though marginal, suggests that the recent decade (D IV) experienced meteorological stressed conditions with respect to the baseline and as compared to D III. Typically, evapotranspiration is limited in such water stress conditions, and is therefore expected to decline in D IV. However, a decreased precipitation not resulting in a decreased ET for D IV suggests that the rapid warming observed in the area somewhat played a key role in the increased actual evapotranspiration. This partly explains why streamflow declined disproportionately in the recent decade.

4.3.2. Decadal Changes in Hydrometeorological Variables

At a glance, the observed increase in precipitation in D III by 9.29% resulting in an increase in flow by 1.46%, and a decrease in precipitation in D IV by 1.83% leading to a decrease in streamflow by 13.22% may lead one to conclude that precipitation is the main climatic factor for streamflow changes in the basin. However, recourse to the elasticity of climate variables to streamflow in the basin will lead to a different conclusion. In

Table 7. Annual streamflow elasticities and contribution of precipitation and actual ET to streamflow changes for respective baseline periods.

		Elasticity (ε)
Period	Contribution to Q Change	Precip	Evapo	Equation
I: 2002~2011
		0.96	−2.88
	relative to BLP I	100%	~	$9.29{\epsilon }_p+2.58{\epsilon }_E=1.46$
		1.35	−0.74
	relative to BLP II	100%	~	$18.13{\epsilon }_p+6.25{\epsilon }_E=19.94$
		1.19	−1.44	~
	relative to BLP III	55%	−30%	$13.54{\epsilon }_p+4.38{\epsilon }_E=9.93$
II: 2012~2021
		0.96	−2.88
	relative to BLP I	13.28%	86.62%	$1.83{\epsilon }_p-3.98{\epsilon }_E=13.22$
		1.35	−0.74
	relative to BLP II	100	~	$6.11{\epsilon }_p+7.71{\epsilon }_E=2.59$
		1.19	−1.44
	relative to BLP III	~	100%	$1.99{\epsilon }_p+5.81{\epsilon }_E=-5.97$

, we find that elasticities of precipitation and actual ET are 0.96 and −2.88, respectively. This suggests that a 10% increase in precipitation results in a 9.6% increase in streamflow, while a 10% increase in actual ET leads to a 28.8% decrease in streamflow. Thus generally, streamflow is less sensitive to precipitation for the reference period. These elasticities also explain the relatively modest (1.46%) increase in streamflow for the 9.29% increase in precipitation in D III, and likewise the 13.22% decrease in flow for an only 1.83% decrease in precipitation in D IV with respect to the 1982–1991 baseline period. Thus, it can be concluded from the above results that actual evapotranspiration is the main climatic factor responsible for streamflow dynamics in the Rockaway sub-basin for this reference period. Although this conclusion holds true, in respect of the actual contribution to the observed streamflow changes for D III, precipitation was entirely responsible with 100% contribution. This means that the amount of precipitation was more than sufficient to satisfy evaporative demands, with the left over going into streamflow generation. For D IV, the impact of actual ET sensitivity to streamflow was largely felt. Such that, while actual ET contributed to approximately 87% of streamflow, precipitation contributed only 13%.

4.3.3. Hydroclimatic Response to Changes Relative to 1992–2001 Baseline, BLP II

Relative to the 1992–2001 reference period, all temperature variables in the basin saw increases consistent with global trends [44,95], with minimum temperature showing the largest increase in D IV (2.25 °C/decade). Unsurprisingly, actual ET followed along with temperature, with increases of 6.25% and 7.71% for D III and IV, respectively. However, these increases were overwhelmed by the respective 18.13% and 6.11% increase in precipitation, leading to a rise in streamflow by 19.94% and 2.59% for D III and IV, respectively (Table 6). Here D III, having precipitation increase by 18.13% can be considered as a meteorologically wet decade compared to D IV and relative to the baseline period. While water limits ET values in dry conditions, energy limits ET values in wet conditions [96,97]. Similarly, actual evapotranspiration in D III was largely energy-limited, leading to a considerable streamflow generation by 19.94%.

In terms of sensitivity to streamflow, Table 7 shows that streamflow elasticity to precipitation is 1.35, indicating that a 10% rise in mean annual precipitation results in a 13.5% increase in streamflow. On the other hand, streamflow elasticity to actual ET is −0.74, which suggests a 7.4% decline in streamflow for a 10% increase in actual ET. This indicates that streamflow is more sensitive to precipitation than actual ET for this reference period and that evaporative demand was overcome by the relative increases in precipitation for D III and IV. As such, the observed increases in streamflow for both decades can be entirely attributed to precipitation as revealed in Table 7.

4.3.4. Hydroclimatic Response to Changes Relative to 1982–2001 Baseline, BLP III

Minimum, mean, and maximum temperatures showed increases in both decades (D III and D IV) relative to the baseline (Table 6). With the recent decade being the warmest in the basin, minimum temperature is observed to have the highest increasing rate (2.37 °C/decade), followed by mean temperature (1.23 °C/decade), and then maximum temperature (0.09 °C/decade). Similarly, precipitation was observed to increase in both decades, with the largest increase in D III, having 170 mm (13.54%) more than in the baseline period. Consistent with temperature, actual evapotranspiration was observed to increase in both decades compared to the baseline. Although the results in Table 6 show that both precipitation and evapotranspiration are increasing for all decades relative to the baseline, they induced varied signal and strength in streamflow. Whereas precipitation was the dominant contributor to streamflow in D III (55%) leading to an increased flow, actual ET entirely overwhelmed precipitation in D IV causing a decrease in flow. The reason for this is that in D III, water was sufficient, and energy becomes the more important control of evapotranspiration whereas in D IV, water was limited, and evapotranspiration, largely driven by energy, went into further decreasing streamflow (−5.97%).

This is confirmed by sensitivity results in Table 7. It reveals that, elasticity of streamflow in relation to precipitation (${\epsilon }_p$) and actual ET (${\epsilon }_{ET}$) for the Rockaway sub-basin are 1.19 and −1.44, respectively. This indicates that 10% increase in precipitation results in 11.9% increase in streamflow while 10% increase in actual ET would lead to 14.4% decrease in streamflow. Thus, annual streamflow was generally more sensitive to the change in actual ET than the change in precipitation, although in D III, precipitation contributed 55% to streamflow changes whereas evapotranspiration contributed 30%.

5. Summary and Conclusions

At a time when the climate change problem has evolved into a crisis, the piece-meal approach to carrying out hydrological impact analysis at a single study location may no longer suffice. At best, a comprehensive study that combines the detection and analysis of trends along with hydrological modelling study will provide important foundations for understanding the hydroclimatic patterns in an area and the driving mechanisms behind these trends in the wake of a changing climate. In this study, we used long-term meteorological and hydrologic observations to identify trends in hydroclimatic indicator variables in the PRB. We also modelled streamflow and groundwater elevation using the Rockaway sub-basin as a case study to understand the impacts of recent climate changes to streamflow in the study basin. Recognizing that hydroclimatic variables, by their nature, are nonmonotonic, we employed the wavelet transform—an advanced trend analysis tool—as against the frequently used MK trend test, to detect and identify patterns in hydroclimatic variables in the PRB. Rather than using the MK trend test directly on the raw data whose results tend to be largely statistically insignificant in respect of precipitation data, for instance, the wavelet transform approach was applied to identify the hidden monotonic trends in the characteristically noisy hydroclimatic time-series. For the hydrological impact assessments, the physically based distributed MIKE SHE hydrologic model provided the platform to successfully simulate the hydrologic conditions of the Rockaway catchment. Based on the model’s water balance outputs, the impacts of recent climate were assessed from changes in naturalized or baseline periods against impacted periods. Further analysis was carried out using climate elasticities to determine the sensitivities and contributions of climatic variables to streamflow changes in the sub-basin to three different baseline conditions. By this, we demonstrated that the time perspective or baseline condition used to assess climate change impacts can also substantially influence results.

Major sources of uncertainty in this study may be that which pertains to hydrological modelling such as input, output, structural, and parametric uncertainties [98,99]. Because streamflow observations used in calibrating and validating the MIKE SHE hydrological model was based on reconstructed data, it is likely that errors emanating from the methods and data used in estimating daily reconstructed streamflow for the Rockaway catchment (refer to [47]) may be propagated in this study. Howbeit, conscious effort was made in minimizing uncertainties in our analyses first by the use of multiple objective function (i.e., observed streamflow and groundwater-level data) that allowed for simultaneous optimization of model parameters. The model’s ability to reasonably simulate both surface and subsurface flows as evidenced in the satisfactory performance criterion give credence to the findings in the study. In addition, one uncertainty that has almost been universally overlooked in climate impact studies is the choice of baseline condition. In our study, we assessed the hydrologic response to changes in climate using three different baseline climates against two recent future periods (i.e., 2002–2011, 2012–2021). We found that the choice of baseline condition constitutes an important source of uncertainty in the sensitivities of streamflow to precipitation and evapotranspiration changes and should routinely be considered in any climate impact assessment. Against this background, we present the key findings from our results below:

(1)

Over the period 1979–2021, minimum, mean, and maximum temperatures showed significantly upward trend in all studied subcatchments of the PRB with minimum temperature having the highest rate of change at 0.059 °Cyr⁻¹ in the Rockaway sub-basin. In contrast, maximum temperatures experienced the slowest rate of change at 0.0034 °Cyr⁻¹. Across the PRB, the rate of change in mean temperature ranges from 0.025–0.035 °Cyr⁻¹.

(2)

Overall, precipitation showed a significant increasing signal in all analyzed sub-basins with the fastest rate of 0.72 mm/yr in the Ringwood catchment and the slowest rate at 0.13 mm/yr in the Rockaway catchment. This observed long-term increasing trend in precipitation and temperature in the PRB is indicative of a changing climate, consistent with the dominant trends in the broader Northeastern region. Spatially, trends in both precipitation intensity (R10) and consecutive dry days (CDD) were observed to decrease in the uppermost portion of the PRB at the Ringwood catchment but increases towards the south in the Rockaway and Upper Passaic sub-basins. This pattern is also dominant in the wider Northeast, and provide further evidence of the connection between extreme weather events and climate change.

(3)

In two out of the three analyzed sub-basins, streamflow displayed significantly downward trends with an increasing trend in the Upper Passaic subcatchment. This is in spite of the increasing trends in both precipitation and temperature in all the three subcatchments. Although it is well established that precipitation amounts and intensity directly affect streamflow [100], the present results rather show that an increase in precipitation does not always lead to an increase in streamflow. From a hydrological modeling standpoint, attempt was therefore made to examine the causes of streamflow in the PRB using the Rockaway subcatchment as a case study.

(4)

Decadal changes in climate revealed that the recent decade (2012–2021) was both the warmest and driest period relative to all baseline periods, and compared with the 2001–2011 decade. It showed a mean temperature increase ranging from 1.16 °C in BLP II and 1.29 °C in BLP I. Being the driest period, the recent decade also showed precipitation changing from −1.83% to 6.11% relative to the 1982–1991 and 1992–2001 baselines, respectively. In contrast, the wettest decade was 2002–2011 relative to all baseline periods with precipitation increase ranging from 9.29% in the 1982–1991 baseline to 18.13% in the 1992–2001 baseline.

(5)

Relative to the overall baseline period (BLP III), the warmest and the driest decade (2012–2021), having a mean temperature increase of 1.23 °C induced an actual evapotranspiration increase of 5.81% and a marginal precipitation increase of 1.99%, resulting in a 5.97% decrease in streamflow. Similarly, the wettest decade (2002–2011), with mean temperature increase of 0.75 °C relative to the overall baseline period (BPL III), induced an actual evapotranspiration increase of 4.38% and a precipitation of 13.54%, which resulted in a streamflow increase of 9.93%.

(6)

Across the three baseline periods, we found that precipitation elasticity to streamflow ranged from 0.96 to 1.35 suggesting that a 10% rise in precipitation will result in between 9.6% to 13.5% increase in streamflow in the study basin. Similarly, evapotranspiration elasticity to streamflow ranged from −2.88 to −0.74 indicating that, a 10% increase in actual ET will lead to between 28.8% to 7.4% decrease in streamflow. The relatively pronounced negative ET elasticity value also reflects the effect of warming climate in the basin. Generally, as temperature increases, ET increases and streamflow decreases. With streamflow showing high sensitivity to actual ET increases more than precipitation, it is safe to conclude that, to a large extent, actual evapotranspiration is more important in the flow dynamics of the PRB in the wake of a warming climate.

(7)

The general observation therefore is that in decades where water is available, energy limits actual evapotranspiration which makes streamflow more sensitive to precipitation increase. However, in meteorologically stressed or dry decades, water limits actual ET thereby making streamflow more sensitive to increases in actual evapotranspiration.

The application of discrete wavelet-transform analysis and process-based hydrological modeling in this study adequately captured the hydroclimatic signatures as well as hydrologic response to climate change in the PRB. A broader study in the future that incorporates how hydrologic sensitivities vary spatially across the PRB will help in further minimizing the uncertainties in climate impact assessments for the basin.