*Article* **Simulating Rainfall Interception by Caatinga Vegetation Using the Gash Model Parametrized on Daily and Seasonal Bases**

**Daniela C. Lopes 1,\*, Antonio José Steidle Neto 1, Thieres G. F. Silva 2, Luciana S. B. Souza 2, Sérgio Zolnier <sup>3</sup> and Carlos A. A. Souza <sup>2</sup>**


**Abstract:** Rainfall partitioning by trees is an important hydrological process in the contexts of water resource management and climate change. It becomes even more complex where vegetation is sparse and in vulnerable natural systems, such as the Caatinga domain. Rainfall interception modelling allows extrapolating experimental results both in time and space, helping to better understand this hydrological process and contributing as a prediction tool for forest managers. In this work, the Gash model was applied in two ways of parameterization. One was the parameterization on a daily basis and another on a seasonal basis. They were validated, improving the description of rainfall partitioning by tree species of Caatinga dry tropical forest already reported in the scientific literature and allowing a detailed evaluation of the influence of rainfall depth and event intensity on rainfall partitioning associated with these species. Very small (0.0–5.0 mm) and low-intensity (0–2.5 mm h<sup>−</sup>1) events were significantly more frequent during the dry season. Both model approaches resulted in good predictions, with absence of constant and systematic errors during simulations. The sparse Gash model parametrized on a daily basis performed slightly better, reaching maximum cumulative mean error of 9.8%, while, for the seasonal parametrization, this value was 11.5%. Seasonal model predictions were also the most sensitive to canopy and climatic parameters.

**Keywords:** rainfall partitioning; dry tropical forest; gash model; interception modelling

#### **1. Introduction**

Water availability is limited in arid and semiarid regions, with rainfall interception playing an important role on site and catchment water balances, as well as in the context of climate change [1]. Rainfall partitioning by trees is an intricate process, mainly affected by canopy and weather factors, such as the characteristics of rainfall events, becoming even more complex where vegetation is sparse [2]. Thus, rainfall interception modelling appears as an important tool for extrapolating experimental results both in time and space, helping to better understand this hydrological process, as well as to implement effective water resource management and land use planning.

Many mathematical models have been developed, validated and successfully applied to simulate rainfall partitioning in different forest types, including coniferous and hardwood stands [3–5], rainforests [6,7], deciduous and sparse canopies [8–10], mixed stands [11] and crops [12]. However, there are few studies about simulating or evaluating the rainfall interception in Caatinga vegetation [1,13–16]. This domain corresponds to an area of tropical dry forest with deciduous tree-shrubs, which covers close to one million km2 in the Northeast of Brazil, occupying around 50% of this region [17]. Caatinga is

**Citation:** Lopes, D.C.; Steidle Neto, A.J.; Silva, T.G.F.; Souza, L.S.B.; Zolnier, S.; Souza, C.A.A. Simulating Rainfall Interception by Caatinga Vegetation Using the Gash Model Parametrized on Daily and Seasonal Bases. *Water* **2021**, *13*, 2494. https:// doi.org/10.3390/w13182494

Academic Editors: Tamara Tokarczyk and Andrzej Walega

Received: 16 August 2021 Accepted: 9 September 2021 Published: 11 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

a fragile ecosystem due to the scarce water resources and the anthropogenic pressures, mainly the intensive exploitation of the region by agriculture and livestock [18]. There is a high temporal and spatial variability of the rainfall regime in this domain, both considering annual and individual events. Specifically, rainfall that occurs in the semi-arid Northeast of Brazil is concentrated over a short period, with the dry season lasting from five to nine months and resulting in uncertainties about the water regime [1]. Caatinga ecosystem is highly dynamic and its vegetation responds quickly to climatic conditions due to morphological and physiological adaptations to aridity by many species of plants. The Caatinga species comprise a whole range of deciduousness, including plants that retain their leaves throughout all the year and other that are leafless during seven months each year [15,17]. The main factor that controls the structure and distribution of vegetation is the precipitation, but photoperiod and nutrients also affect the Caatinga species [18].

More studies about rainfall interception in the Caatinga domain are important to increase the spatial and temporal accuracy in rainfall partitioning simulations and to better understand this process in dry tropical forests. Additionally, these studies may also benefit watershed and forest managers of other similar arid and semiarid ecosystems, since dry tropical forests are recognized as one of the world's major biomes and are found in a wide area extending from the Amazon basin in South America towards northern Mexico and the Caribbean [19].

The analytical Gash models [9,20] are most often used when predicting rainfall interception due to their ease of use, low parameter requirement and low programming complexity [21]. These models are capable of estimating rainfall partitioning by using a series of parameters based on canopy structure, evaporation rate and rainfall regime [22]. The original analytical Gash model [20] represents rainfall input as series of discrete storms, each comprising a wetting up period, a saturation period and a drying out period [7]. The sparse version [9] encompasses the case of forest stands with significant open spaces between tree canopies, also introducing some minor corrections [21]. The main difference between these two versions is that the sparse model is based on evaporation and canopy storage per unit area of canopy cover rather than per unit of ground area. This overcame a limitation in the description of sparse forests by the original model, which can prevent the simulated canopy from wetting up [9,11].

Both original and sparse Gash models [9,20] are typically applied using mean annual or seasonal rainfall intensity and evaporation rates, which are considered as constant parameters in all events during the simulated period [2,8,13]. The same occurs with the canopy storage capacity, canopy cover fraction and threshold value required to saturate the canopy [5,23,24]. The sparse Gash model was already parametrized on a daily basis, considering a linear relationship between leaf area index and canopy storage capacity during the plant cycle [12,22]. These model adaptations, based on estimates of parameters for individual storm events, were also compared to other methodologies [3], reinforcing that the daily changes, observed in canopy structures, especially for deciduous vegetation, tend to reduce systemic simulation inaccuracies.

Although a number of works have been focused on applying and evaluating the Gash model with different parametrizations and for distinct forest types, its parametrizations on daily and seasonal bases were not studied for the Caatinga domain. These procedures tend to better represent the effects of changes in canopy cover on the rainfall interception process, mainly in this deciduous ecosystem, where canopy structure often changes gradually, but relatively rapidly. Furthermore, the adjustments and modifications required to perform such simulations allow a more detailed and accurate evaluation of the influence of rainfall depth and event intensity on rainfall partitioning. Therefore, the objective of this study was to parametrize the sparse Gash model on daily and seasonal bases, validating these approaches for simulating rainfall interception from five Caatinga species (*Spondias tuberosa*, *Commiphora leptophloeos*, *Cnidoscolus quercifolius*, *Aspidosperma pyrifolium and Cenostigma pyramidale*), improving the description of this hydrological process already reported in the scientific literature for dry tropical forests by enhancing the temporal and spatial accuracy of the estimates performed by the model. Additionally, the influence of rainfall depth and event intensity on rainfall partitioning associated with these species was evaluated.

#### **2. Materials and Methods**

#### *2.1. Meteorological and Rainfall Measurements*

This study was conducted within a private property with Caatinga vegetation area of 81,000 m2 and density of 930 trees ha−1, located in the Floresta municipality, Pernambuco State, Brazil (08◦18 31" S, 38◦31 37" W, 378 m a.s.l.). Vegetation in the experimental plot is mainly composed by the native species *Spondias tuberosa*, *Commiphora leptophloeos*, *Cnidoscolus quercifolius*, *Aspidosperma pyrifolium* and *Cenostigma pyramidale* (Figure 1), which are randomly distributed over the study site and are representative of the Caatinga domain [17,25–27]. Table 1 presents the main characteristics of these tree species [19], which are common to other species found in the Caatinga domain [14,16].

**Figure 1.** Study area, monitored trees and micrometeorological tower.


**Table 1.** Main characteristics of the tree species used in the rainfall interception measurement: number of individuals (N), average diameter at breast height (DBH), average number of stems (NB), average tree height (H) and average tree crown projected area (CPA).

The climate of the region according to the Köppen classification is BSh, corresponding to a tropical semiarid hot type [28], and the Thornthwaite Aridity Index is 0.48 [29], confirming the region as semiarid. Average annual precipitation, wind speed and net solar radiation are 489.3 mm, 2.3 m s−<sup>1</sup> and 22.2 MJ m−<sup>2</sup> day−1, respectively. Mean annual air temperature is around 26.1 ◦C, with mean monthly temperatures ranging from 23.3 ◦C in July to 28.3 ◦C in November. Mean annual relative humidity is 61.9%, with mean monthly values between 50.3% in October and 70.4% in April.

Micrometeorological measurements were performed by electronic sensors installed on a galvanized iron tower at 8 m above the ground (Figure 1). Data were registered and stored in a datalogger (CR10X, Campbell Scientific, Logan, UT, USA) and measurements were conducted continuously from 1 March 2016 to 30 September 2017. The wind speed and direction were measured by an anemometer (03002 Wind Sentry, R. M. Young Company, Traverse, MI, USA). A quantum sensor (SQ-321, Apogee Instruments, Logan, UT, USA) measured the photosynthetically active radiation (PAR), while radiation balance was measured by a net radiometer (NR-Lite, Kipp & Zonen, Delft, The Netherlands) and the global solar radiation was obtained by a pyranometer (SP-230, Apogee Instruments, Logan, UT, USA). An automatic rain gauge (CS700-L, Hydrological Services Pty, Sydney, Australia) registered the gross rainfall.

The photosynthetically active radiation transmitted through the canopies was measured at two below-canopy positions, previously defined in representative trees of the predominant Caatinga species in the study area, by hand-moving two linear quantum sensors (SQ-321, Apogee Instruments, Logan, UT, USA) from one tree to another. Additionally, three aspirated psychrometers, made of T-type thermocouples (copper-constantan), were used for obtaining the dry and wet bulb temperatures at 0.5, 1.5 and 2.5 m above the mean canopy level. Two soil heat flux plates (HFT3-REBS, Hukseflux, Delft, The Netherlands) were also installed in the top soil layer at a depth of approximately 0.05 m.

A ceptometer (LP-80, Decagon Devices, Pullman, DC, USA) was used to measure the fractional interception of photosynthetically active radiation. The incident radiation measurements were performed in open areas without physical obstacles, not including cloudy days or at dusk. The transmitted radiation measurements occurred under the tree canopies. One incident and four transmitted radiation measurements in different directions (north, east, south and west) were executed for each sample of predominant Caatinga species, during 14 campaigns with 135 readings each. Based on the fractional interceptions of photosynthetically active radiation, the integrated microprocessor of the ceptometer estimated the leaf area index based on a simplified version of the Norman– Jarvis radiation transmission and scattering model [30,31]. Polynomial equations were then fitted to measure fractional interception of photosynthetically active radiations obtained by the ceptometer and registered in the datalogger for estimating daily leaf area indices for each studied species.

Throughfall was measured by 15 manual collection gauges, placed randomly underneath the vegetation canopy, comprising three gauges per predominant Caatinga species. Measurements were performed after each rainfall event and the gauges were installed at 1.0 m above ground level, presenting orifices of 0.07 m2. Gauges were installed at a halfway distance between canopy edge projection and stem in order to minimize the effects of spatial variability on the magnitude of average throughfall [13]. The area under each tree was divided by three diagonals considering the crown projected limits, totaling six sampling points equally spaced at angles of 60◦. The gauges were distributed in three sampling points, being representative of a 120◦ circular sector bisected by the gauge longitudinal axis and centered on the tree position. Each gauge was relocated after every three rainfall events to a new position correspondent to the empty sampling point located beside it and following the area clockwise. This procedure minimizes errors originating from spatial variability and improves long-term sampling [32,33]. Furthermore, it allows to derive reliable mean throughfall per tree, even with a limited number of gauges [2,10,23,34].

Stemflow was measured after each rainfall event by installing twelve zinc gutters of 0.15 m in height, attached to the tree stems at 1.3 m above ground level and connected to individual plastic containers. A hose was fed into each plastic container and measurements were performed with a graduated test tube, with the purpose of reducing evaporation. Due to the tortuous trunks and rough bark of *S. tuberosa* and *C. leptophloeos*, the stemflow monitoring was restricted to the other three species, which presented a projected crown radius greater than 0.2 m and were sub-divided into two classes of diameter at breast height [35], that were 0.05 m ≤ DBH < 0.10 m and 0.10 m ≤ DBH ≤ 0.20 m. This sampling plan assured reliable stemflow measurements, since the two unmonitored trees have canopy structures similar to other species with low stemflow [36]. *Spondias tuberosa* presents low inclined branches with many flow path obstructions that create drip points, enhancing throughfall production, while *Commiphora leptophloeos* has horizontal leaves and only one stem.

The effective rainfall interception was obtained for each rain event by subtracting the gross rainfall by the sum of measured throughfall and stemflow [34,37]. Each rainfall event was defined as a period when cumulative gross rainfall exceeded 0.2 mm, provided that there was a minimum of 6 h without rainfall between events [2,7,24,38]. A consistency analysis was performed on the rainfall and micrometeorological data with electronic spreadsheet functions to remove all inconsistent values and outliers. Visual analysis of graphs relating the variables to time complemented the data evaluation.

#### *2.2. Sparse Gash Model Parametrized on Daily and Seasonal Bases*

A daily parametrization for the sparse Gash model was proposed in this study for simulating the rainfall interception of Caatinga species [3,12,22] and was compared with a seasonal parametrization (Table 2), which was based on mean constant parameters for rainy and dry seasons.

In both approaches, the net rainfall interception was estimated as [9]

$$\text{IN} = \text{IC} + \text{IW} + \text{IS} + \text{IA}, \tag{1}$$

where IN is the net rainfall interception (mm), IC is the interception insufficient to saturate the canopy (mm), IW is the rainfall interception during canopy wetting (mm), IS is the rainfall interception during saturated canopy conditions (mm) and IA is the evaporation after rain ceased (mm).

The original and sparse Gash models also include a formulation for stemflow and evaporation of water stored on wetted trunks. However, the cumulative stemflow accounted for only 0.7% of total rainfall for the five studied Caatinga species and was considered negligible during simulations in this work. Thus, an equivalent interception was calculated, corresponding to the difference between gross rainfall and throughfall [7,8,34].


**Table 2.** Equations describing the components of rainfall interception in the seasonal and daily parametrizations of the sparse Gash model proposed in this paper.

M, number of storms insufficient to saturate the canopy (dimensionless); i, mean value for a rainfall event (dimensionless); c, canopy cover fraction (dimensionless); y, mean value for rainy or dry season (dimensionless); PG, gross rainfall (mm); n, number of storms sufficient to saturate the canopy (dimensionless); PS, threshold value required to saturate the canopy (mm); Sc, canopy storage capacity per unit area of cover (mm); Ec, evaporation rate from wet canopy per unit area of cover (mm h−1); R, rainfall rate or rainfall intensity (mm h−1).

#### *2.3. Estimation of Meteorological and Canopy Parameters*

The threshold value required to saturate the canopy (PS, mm) was obtained on a seasonal or daily basis, depending on the model approach, as [9,12,22]

$$P\_{\rm Sy} = -\left(\frac{\rm R\_{y}S\_{\rm Cy}}{E\_{\rm Cy}}\right) \ln\left(1 - \frac{E\_{\rm Cy}}{R\_{\rm y}}\right),\tag{2}$$

$$P\_{\rm Si} = -\left(\frac{\rm R\_i S\_{\rm Ci}}{E\_{\rm Ci}}\right) \ln\left(1 - \frac{E\_{\rm Ci}}{R\_i}\right),\tag{3}$$

When using the parametrization on a daily basis, the rainfall rate (Ri, mm h<sup>−</sup>1) was the average rainfall intensity during all hours in each storm event [22]. For the seasonal sparse Gash model, mean rainfall rates (Ry, mm h−1) were calculated separately [24], for rainy (December–May) and dry (June–November) seasons and then applied in a generalized form to all individual rainfall events. Rainy and dry periods were determined according to the rainfall pattern observed in the studied region, as well as the phenological and leaf area index data of the five studied Caatinga species [1,16,19].

The evaporation rate from wet canopy (Em, mm h−1), which represents the evaporation from the canopy during the storms, was also calculated for each storm event for the parametrization on a daily basis, while mean values were obtained for rainy and dry seasons when using the seasonal approach. This parameter was estimated hourly based on the Penman–Monteith equation [39], excluding the non-storm periods, with the canopy resistance set to zero [40] and using the momentum method for estimating the aerodynamic resistance [5,39]. The estimated evaporation rates from wet canopy were then divided by the canopy cover fractions before being applied in the models (Ec, mm h<sup>−</sup>1), adjusting the original values for a complete canopy in proportion to the canopy cover [9].

The canopy cover fraction (c, dimensionless) described the vegetation density. Daily canopy cover fractions were calculated according to the Beer–Lambert equation [41]:

$$\mathbf{c} = \mathbf{1} - \mathbf{P} \mathbf{B} / \mathbf{P} \mathbf{A} = \mathbf{1} - \mathbf{e}^{-\mathbf{k} \cdot \mathbf{L}},\tag{4}$$

where PA is the incoming photosynthetically active radiation on canopy (μmol m−<sup>2</sup> s−1), PB is the transmitted photosynthetically active radiation through the canopy (μmol m−<sup>2</sup> s<sup>−</sup>1), k is the extinction coefficient (dimensionless) and L is the leaf area index (m2 m<sup>−</sup>2).

The Beer–Lambert model describes the radiation transmittance through crop canopies as an exponential-type attenuation process, which can be also associated with the fractional photosynthetically active radiation, as well as with the leaf area index [12]. Daily canopy cover fractions were directly applied in the sparse Gash model parametrized on a daily basis, while average values for rainy and dry seasons were obtained when performing simulations with the seasonal parametrization.

The canopy storage capacity (S, mm) corresponded to the amount of water remaining in the canopy after rainfall and throughfall cease, considering evaporation equals to zero [42,43]. This parameter depends on the intensity and duration of the storm, as well as the spatial and temporal variability of trees [44]. In this study, S was assumed to have a linear relationship with the leaf area index [12]. Furthermore, the canopy storage capacity was adjusted per unit area of cover (Sc, mm), by dividing the original S value by the canopy cover fraction before applying it in the models [9].

The mean method was used for estimating a specific S value [3,45], representing the depth of water retained by leaves of each studied species. For this, scatter plots of measured rainfall interception versus gross rainfall were plotted for a number of rain events large enough to saturate the canopy of each Caatinga species and the specific canopy storage capacities were derived from the intercepts of the regression lines fitted to these data. That is, two regression lines were created, relating rainfall interception to gross rainfall for storms that are either insufficient or sufficient to saturate the canopy. The slope of each regression line was determined by an iterative least square fitting procedure. The difference between gross rainfall and throughfall at the intersection point of these two regression lines provided the estimate of S. The use of rainfall interception when plotting the regression lines yields the least stochastic errors, mainly when rainfall outside the canopy is measured without observational scatter and rainfall inside the canopy is observed with scatter [13,45]. For the simulations with daily parametrization, specific S values were multiplied by the daily leaf area index of each species, resulting in different daily S values. When applying the seasonal parametrization, daily S values were averaged considering rainy and dry periods.

#### *2.4. Validation Analysis*

The predictive capacity of the adjusted model was evaluated by the statistical parameters cumulative mean relative error, mean bias error, index of agreement and Nash–Sutcliffe efficiency [23,46,47]:

$$\text{CMRE} = 100 \frac{|\mathbf{C}\_{\text{I}} - \mathbf{C}\_{\text{S}}|}{\mathbf{C}\_{\text{S}}},\tag{5}$$

$$\text{MBE} = \frac{\sum \left( \mathbf{P}\_{\mathbf{j}} - \mathbf{O}\_{\mathbf{j}} \right)}{\mathbf{w}},\tag{6}$$

$$\mathbf{d} = 1 - \frac{\boldsymbol{\Sigma} \left( \mathbf{P}\_{\mathrm{j}} - \mathbf{O}\_{\mathrm{j}} \right)^{2}}{\boldsymbol{\Sigma} \left( \left| \mathbf{P}\_{\mathrm{j}} - \mathbf{O}\_{\mathrm{m}} \right| + \left| \mathbf{O}\_{\mathrm{j}} - \mathbf{O}\_{\mathrm{m}} \right| \right)^{2}} \tag{7}$$

$$\mathcal{E} = 1 - \frac{\sum \left( \mathcal{P}\_{\text{j}} - \mathcal{O}\_{\text{j}} \right)^{2}}{\sum \left( \mathcal{O}\_{\text{j}} - \mathcal{O}\_{\text{m}} \right)^{2}} \tag{8}$$

where CMRE is the cumulative mean relative error (%), CI is the real cumulative rainfall interception (mm), CS is the simulated cumulative rainfall interception (mm), MBE is the mean bias error (mm), Oj is the measured rainfall interception (mm), Pj is the predicted rainfall interception (mm), w is the number of testing data (dimensionless), d is the index of agreement (dimensionless), Om is the average experimental rainfall interception (mm) and E is the Nash–Sutcliffe efficiency (dimensionless).

Values of MBE close to zero indicate that the model is useful for prediction, with negative and positive values suggesting underestimates and overestimates, respectively [46]. This indicator is related with the unit in which the evaluated property is expressed, as well

as with the dataset range of values, which, in this study, represents the rainfall interception, in mm.

The CMRE, d and E values are standardized measures in which cross-comparisons for a variety of models, regardless of units, are possible. That is, these indicators are not measures of correlation or association in the formal sense, but rather measures of the degree to which the predictions obtained from a model are error-free [47,48]. E and d vary from 0 to 1, with the maximum value representing a perfect agreement between observed and predicted data. CMRE ranges between 0 and 100%. Based on approximately 111 scientific research studies about rainfall interception modelling [21], the cumulative mean error was classified in five qualitative groups: bad (>30%), applicable (10–30%), good (5–10%), very good (1–5%) and extremely good (<1%).

Additionally, validation graphs of the measured rainfall interceptions against the predicted ones were plotted. Aiming at verifying the model accuracy, the t-test was applied to the intercept (b) of each linear regression to check whether it was significantly different from 0 and to the line angular coefficient (a) to confirm whether it was significantly different from 1, at the level of 1% probability.

#### *2.5. Statistical and Sensitivity Analyses*

Sensitivity analyses were performed to identify the relative importance of canopy and climatic parameters (S, c, Em and R) in both daily and seasonal parametrizations of sparse Gash model. For this, the values of each parameter were increased and decreased by up to 50% of their original values and the simulated results were compared to measured data [5,11,24].

Measured rainfall interceptions and estimated model parameters were also submitted to variance analysis and averages were compared by the F and Scott–Knott tests (*p* < 0.05) through the SISVAR software (Federal University of Lavras-UFLA, Lavras, Minas Gerais, Brazil) [49]. Statistical differences between Caatinga species, storm classes and simulation periods were evaluated.

#### **3. Results**

#### *3.1. Rainfall Partitioning*

The total measured rainfall between 1 March 2016 and 30 September 2017 was 429.5 mm, generated by 66 discrete rainfall events. From these, 343.7 mm (80.0%) occurred during the rainy season, while 85.8 mm (20.0%) occurred during the dry season. The average, maximum and minimum event-based rainfall amounts were 6.5 (±9.3), 36.4 and 0.2 mm, respectively. Rainfall intensities varied from 1.2 to 19.2 mm h<sup>−</sup>1, with an average of 3.2 (±2.9) mm h<sup>−</sup>1.

Statistical analyses indicated that the frequency distributions of the event size (Table 3) and intensity (Figure 2) among annual analysis, rainy and dry seasons did not differ significantly. The very small events (0.0–5.0 mm) were significantly more frequent than the other four classes. However, as shown in Table 3, they contributed with the lowest percentages to total gross rainfall during the rainy season and for the annual analysis. When evaluating the dry season, events from 20.1 to 40.0 mm were not observed and small events (5.1–10.0 mm) were responsible for the lowest percentages of total gross rainfall. The highest percentages of total gross rainfall were verified for the very large events (30.1– 40.0 mm) during the rainy season and for middle events (10.1–20.0 mm) during dry and annual periods.


**Table 3.** Number of events (NE), cumulative gross rainfall (CGR) and percentage of gross rainfall (PGR) in different rainfall classes during the study period.

**Figure 2.** Frequency distributions (bars) and percentages of gross rainfall of the intensity of rainfall events (lines) at the Caatinga domain in Brazil during the measurement period (annual, rainy and dry seasons) from 1 March 2016 to 30 September 2017.

When evaluating rainfall intensity, events from 0 to 2.5 mm h−<sup>1</sup> presented significantly higher frequency of occurrence. Figure 2 shows that these low intensity rainfall events contributed to a higher percentage of gross rainfall only during the dry season, while middle and large intensity rainfall events (2.6–10 mm h<sup>−</sup>1) were responsible for the greatest percentages of gross rainfall both for annual and rainy season conditions.

As shown in Figure 3, the annual cumulative interception values were between 90.5 and 169.7 mm, resulting in average proportions of gross rainfall into interception of 27.9% (±6.5). When considering the rainy and dry seasons, the cumulative interceptions ranged from 43.5 to 114.7 mm and from 35.6 to 55.0 mm, with proportions of gross rainfall to interception of 17.2% (±5.7) and 10.7% (±1.8), respectively. Additionally, the annual average proportions of gross rainfall to throughfall and stemflow were 71.6% (±7.5) and 0.7% (±0.2), respectively.

**Figure 3.** Cumulative measured gross rainfall and interception values from March 2016 to September 2017 at the Caatinga domain in Brazil.

Observing the different storm size classes (Table 4), the proportions of gross rainfall into interception of very small storms (<5 mm) were significantly higher for all studied Caatinga species. Statistical analysis also showed that there was no significant difference between large (20.1–30.0 mm) and very large (30.1–40.0 mm) storms, as well as between small (5.1–10.0 mm) and middle (10.1–20.0 mm) storms. Filtering by significantly equal rainfall classes (<5, 5.1–20.0, 20.1–40 mm), the proportions of gross rainfall to interception tended to decrease as gross rainfall increased. For very large events (30.1–40 mm), *C. quercifolius* presented a significantly lower percentage of rainfall interceptions (Figure 3 and Table 4) and, for the other storm classes, there was also a trend of smaller values when compared with the other studied species.

#### *3.2. Model Parameters*

The leaf area indices, as well as the parameters for the adapted Gash model applied to the five studied Caatinga species, are shown in Table 5. The values of each parameter and leaf area index did not significantly differ among Caatinga species, but *S. tuberosa* presented a larger leaf area index, Sc and c, which also helps to explain the significantly higher rainfall interception observed during the experimental trial. Parameters and leaf area index were statistically equal between seasons and annual analysis, but their numerical differences contributed to improve the simulation results, as discussed below.


**Table 4.** Proportions of gross rainfall portioned into interception (I:GR) considering five rainfall classes.

#### *3.3. Sensitivity Analyses*

The sensitivity analysis of the canopy and climatic parameters to rainfall interception is presented in Figure 4. Variations observed for the five studied species were averaged for each model, since they presented very similar patterns.

**Figure 4.** Sensitivity analysis of the Gash model parametrized on daily and seasonal bases applied to the Caatinga domain in Brazil.

**Table 5.** Measured leaf area indices (L) and estimated canopy storage capacity per unit area of cover (Sc), relative evaporation rate per unit area of cover (Ec/R), evaporation rate from wet canopy (mm h−1), canopy cover fraction (c) and threshold value required to saturate the canopy (PS) for the sparse Gash model parametrized on daily and seasonal bases applied to each studied Caatinga species. Values in brackets represent standard deviations.


Canopy storage capacity (S), canopy cover fraction (c) and evaporation rate from wet canopy (Em) presented positive relationships with interception, whereas mean rainfall rate (R) resulted in a negative relationship with interception. All parameters were sensitive to interception, but S had larger effects on the Caatinga rainfall partitioning, followed by Em, c and R, respectively. Seasonal model predictions were most sensitive to canopy and climatic parameters, which were more pronounced for the S variable.

The sensitivity analysis showed that Em changes could lead from −13.8 (±1.5) to 13.4 (±1.3)% errors in rainfall interception when applying the sparse Gash model parametrized on a daily basis, while, for the seasonal approach, errors could vary from −14.6 (±1.6) to 14.5 (±1.6)%. Considering the R parameter, a −50% change in daily and seasonal parametrizations caused an average 3.2 (±1.5)% difference in the simulated rainfall interception.

For the daily parametrization, a decrease of 50% in c resulted in an average decrease of 14.1% (±0.8) in simulated interception, while an increase of 50% resulted in an average rise of 7.3% (±0.5). When applying the seasonal parametrization, a decrease of 50% in c could lead to an average decrease of 29.7% (±3.4) in simulated interception, while an increase of 50% could lead to an average rise of 17.0% (±0.9).

If the value of S increased by 50%, simulated rainfall interception tended to rise by 25.7% (±1.2) and 59.6 (±3.2) on average by applying the daily and seasonal parametrizations, respectively. On the other hand, the decrease of 50% resulted in average reductions of 28.8% (±1.4) and 43.3% (±2.3) for daily and seasonal parametrizations, respectively.

#### *3.4. Rainfall Interception Simulations*

Figure 5 shows that the daily and seasonal parametrizations performed very similarly, when simulating cumulative rainfall interception, with the daily parametrization performing slightly better, except for *S. tuberosa*. The average annual proportions of gross rainfall into interception were simulated as 27.1% (±5.1), when applying the daily parametrization, and as 26.5% (±5.3), when using the seasonal parametrization. Compared with the

measured data, the differences were 0.8 and 1.4%, respectively. When considering the tree species, simulations for *A. pyrifolium* and *C. pyramidale* resulted in more accurate estimates. These species reached cumulative relative mean errors from 1.23 to 3.69%, considered very good [21]. The cumulative relative mean errors for the other species varied from 5.00 to 8.98% and were considered good [21].

**Figure 5.** Cumulative measured interception values compared with those estimated by the sparse Gash model parametrized on daily and seasonal bases for each Caatinga species in Brazil.

On a per storm analysis (Figure 6), when estimates were evaluated for individual rainfall events, both simulation approaches also performed similarly, with the daily

parametrization presenting slightly less scatter. Prediction errors were higher mainly for small, middle and large storms (5.1–30 mm), with interceptions between 0.2 and 6 mm. Among the tree species, best results were obtained with *S. tuberosa*.

**Figure 6.** Validation graphs of experimental interception values against those predicted by the sparse Gash model parametrized on daily and seasonal bases applied to Caatinga species.

The reliability of the proposed models for Caatinga species was also proven by the statistical indicators (Table 6). MBE varied from −0.20 to 0.15 mm for daily and seasonal parametrizations, indicating predicted values close to measured ones and confirming the model trends of slightly underestimating the *S. tuberosa* interception for both tested approaches, as well as for *A. pyrifolium* with seasonal parametrization and *C. pyramidale* with daily parametrization. The other cases showed trends of slightly overestimating. The d and E averages were 0.94 (±0.04) and 0.75 (±0.80), respectively, when applying seasonal simulations. For the daily parametrization, these indices reached 0.94 (±0.03) and 0.76 (±0.12) on average, respectively. The d and E results proved the high accuracy and the very good agreement between measured and predicted interceptions for all Caatinga species, as well as the slightly better performance of the daily parametrization.

**Table 6.** Summary statistics of interception values predicted by the sparse Gash model parametrized on daily and seasonal bases applied to Caatinga species on storm-based rainfall analysis.


#### **4. Discussion**

#### *4.1. Rainfall Partitioning*

The partitioning pattern observed in this study for the five Caatinga species agreed with other studies of semiarid regions [1,2,16,24], where rainfall is concentrated over a short period. Highly variable rainfall depths per event (0.2–40 mm) were also reported when applying the Gash model to deciduous shrubs in the semiarid Qinghai–Tibet Plateau [22]. Furthermore, mean values of gross rainfall per event (7.2 and 5.1 mm) close to the average value found in this work were verified for semiarid regions of Spain and Iran, respectively [2,23]. On the other hand, average Caatinga rainfall intensities (3.2 mm h−1) were higher than those observed in other semiarid regions. For example, a mean rainfall intensity of 1.7 mm h−<sup>1</sup> was found when modelling the interception in central–western Spain [2], while rainfall intensities equal to 1.8 and 2.3 mm h−<sup>1</sup> were reported for semiarid regions of China and Kenya, respectively [24,50]. The frequency distributions of the event size and intensity were in agreement with other studies about rainfall partitioning in semiarid regions [16,23,24].

Average rainfall interception was significantly higher during the rainy season, which is justified by the largest number of rainfall events and greatest rainfall amounts observed during this period, but also by the reduction in leaf amounts, as well as canopy cover, during the dry season. Rainfall interception of *S. tuberosa* was significantly higher than those verified for the other studied Caatinga species. This difference can likely be related with canopy characteristics, mainly the number of stems, diameter at breast height and tree crown projected area (Table 1), on which further studies are required. Additionally, the leafless periods of *S. tuberosa* corresponds to around 3 months [26], while the other species are leafless from approximately 4 to 7 months [25,27,51] and the emergence of leaves strongly affects the interception process, modifying the redistribution by the tree and the profile of rainwater.

The trend of decreasing the proportions of gross rainfall to interception as gross rainfall increases was also verified in other studies [23]. Additionally, the lower percent rainfall interceptions of *C. quercifolius* can be explained since the peak of its leaf fall lasts around 5 months [25] and this tree produces small leaves during the dry season that only attain their maximum size during the rainy season. For this, it presents a less dense canopy, which tends to facilitate water flow in throughfall [15].

Results obtained in this study confirm the findings of other scientific research studies [1,8,16,23], which suggested that changes in the proportions of gross rainfall to interception are mainly associated with the size of gross rainfall and can be explained since most of the gross rainfall is stored in the canopy during the very small and small rainfall storms. On the other hand, canopies tend to saturate during large and very large rainfall events, increasing throughfall, while the remaining rainfall is stored in the canopies and lost as evaporation during the storm event [9].

#### *4.2. Model Parameters*

Average leaf area indices observed in the Caatinga domain were similar to those verified for species of other semiarid regions, such as those for *R. pseudoacacia* in the Shaanxi province, China [52], as well as for *Q. ilex* and *Q. pyrenaica* in Sardon stands, Spain [2]. However, *S. tuberosa* stood out in this study, maintaining most of its leaves during the dry season, as well as reaching higher individual leaf area indices, when compared with other Caatinga species and native trees of different semiarid regions [22,23]. The average leaf area indices for rainy and dry seasons (Table 5) reflected the seasonal variations of this parameter, which increased linearly from the beginning to the end of the rainy season and decreased linearly during the dry season for all studied species. Indeed, the peak of leaf flush for the Caatinga species tends to coincide with the rainy season, but this process is also affected by the photoperiod [25].

Despite rainfall interception is known to be closely related to the leaf area index [53], other factors influence this process. The increase in canopy density causes leaves to touch, hindering the fully saturation of the entire canopy [12]. Furthermore, wind may reduce canopy storage capacity, as well as branch shape, and leaf inclination and canopy thickness may turn the leaves less wettable.

As expected and as shown in Table 5, both c and Sc followed the pattern observed in the leaf area index for all studied species, with average values of the rainy season larger than those of the dry season. The same behavior was verified for the PS values, since all of these parameters derived from the leaf area index in the parametrizations performed in this study. The c values were similar to other simulations of rainfall interception for deciduous forests in semiarid regions [8,14,22], but Sc and PS were larger. These results are consistent with the larger leaf area indices verified mainly during the rainy season and are probably also associated with the height of the studied Caatinga species (Table 1). The Sc and PS values obtained in this study can be also justified by the distinct forest structures among semiarid regions, since the canopy morphology and physiology interfere in the parametrization, including the leaf amounts and canopy cover during leafed and leafless periods [17]. Additionally, Sc and PS are closely related to Caatinga weather conditions, mainly rainfall distribution and intensity [38]. The estimate of Sc also depended on the specific canopy storage capacity, which was based on the mean method, relating the observed gross rainfall, interception and throughfall that were either insufficient or sufficient to saturate the canopy. On the other hand, the PS also depended on the evaporation rate from wet canopy, which was estimated by the Penman–Monteith equation, considering meteorological (air temperature, net radiation, vapor pressure deficit and wind speed) and vegetation characteristics (crop height).

Canopy coverage is an important parameter in both daily and seasonal parametrizations, since it is a structural parameter, directly related to the free throughfall coefficient and the leaf area index [5]. The free throughfall coefficient, which is assumed to be one minus canopy cover, affects soil water content and nutrient cycling, since it reflects the fraction of rainfall passing through the canopy without contacting the tree surface or removing dry leaves and twigs in the canopy [24]. Indeed, Table 5 shows a trend of species with smaller leaf area indices and greater leafless periods (*C. pyramidale* and *C. quercifolius*) presenting

lower canopy coverages, with consequent larger free throughfall coefficients. The obtained c values can also explain the lower rainfall interceptions verified for *C. quercifolius*, as well as the significantly high interceptions from *S. tuberosa*. These species presented respectively the lowest (0.27) and the largest (0.97) canopy coverages when applying the daily parametrization, with *C. quercifolius* resulting in the greatest range of c values during the study and the smallest canopy coverage during the dry season (Table 5). Additionally, *S. tuberosa* was the species with the largest c values both in dry and rainy seasons. Thus, there is a trend of *S.tuberosa* species protect the Caatinga floor from raindrop splash erosion, also delaying the peaks in storm runoffs, as was verified in semiarid Northeast of China [24] and Brazil [54]. On the other hand, these trees could contribute to the enhancement of soil water scarcity in the Caatinga domain, since the high c values lead to less throughfall reaching the forest floor. However, detailed comprehension about these effects merit further studies.

Considering the daily parametrization, Em was highly changeable during the studied events, varying between 0.19 and 0.94 mm h−1, with an average of 0.40 (±0.14) mm h−1. The Em values observed during rainy and dry seasons reached 0.41 (±0.02) and 0.39 (±0.02) mm h−<sup>1</sup> on average, respectively. Scaled to canopy cover, Ec values ranged from 0.24 to 2.50 mm h<sup>−</sup>1, with an average of 0.67 (±0.36) mm h−<sup>1</sup> for the daily parametrization, while the seasonal approach resulted in Ec values of 0.62 (±0.08) and 0.57 (±0.05) mm h−<sup>1</sup> on average during rainy and dry seasons, respectively (Table 5). Maximum in-storm evaporation rates from 1.83 to 3.98 mm h−<sup>1</sup> were also observed in tropical dry and semiarid regions of Mexico [47], but the average evaporation rates found in this study were larger than those of semiarid regions of Spain, Iran and China [2,23,24]. These results are related with the semiarid climate type of Caatinga, which is mainly characterized by high temperatures and solar radiation, tending to increase evaporation rates when compared with other semiarid and dry tropical stands [1]. The higher evaporation rates and the distinct rainfall intensities affected the Ec/R ratios, which were lower during the rainy season (0.15 ±0.01), reaching average values of 0.26 (±0.03) during the dry season. When applying the daily parametrization, Ec/R presented minimum values between 0.03 and 0.04, while maximum values varied from 0.79 to 0.86 for the five studied species.

The differences observed between seasonal and daily parametrizations, regarding the sensitivity analysis, are expected, since the seasonal model uses two sets of constant parameters, while the daily parametrization considers the daily changes of canopy structure and weather conditions, tending to better represent the associated components and processes of rainfall interception. The sensitivity analysis agreed with other studies [53,55], which found that canopy storage capacity is among the most influential parameters on simulated rainfall interception in deciduous vegetation. Other factors that affect the sensitivity of model predictions are the rainfall and climate characteristics, such as raindrop size distribution, rainfall intensity and wind speed, though those factors are not included in the Gash model [24]. Additionally, canopy basal area and height interfere in the parametrization process, as well as the woody light-blocking elements from the canopy with respect to diameter growth, represented by the wood area index [56].

#### *4.3. Rainfall Interception Simulations*

The good results observed when simulating cumulative rainfall interception may be consequence of using leaf area indices during estimates of c and S, as well as the simulations based on daily or seasonal parameter variations, which allowed the proposed models to describe the rainfall interception patterns better than other approaches [13,14].

When considering the per-storm simulations (Figure 6), very large rain events (30.1–40 mm) were less frequent (Table 3), representing 6.1% of total gross rainfall, and were not observed during the dry season. This probably led to less events with interceptions greater than 6 mm and, consequently, less scatter was verified for these situations. There were more outliers for smaller interceptions, with simulated values moving further from the 1:1 line in both parametrizations and showing that model estimation was less accurate in some individual events. For these cases, the daily parametrization was also subtly better than the seasonal one. These discrepancies certainly affected the statistical indicators (Table 6). However, the scatter patterns agreed with other works [8,23] and, despite the observed discrepancies when comparing measured and simulated individual interceptions, both daily and seasonal parametrizations resulted in a good fit, with all intercept and angular coefficients not significantly differing from 0 and 1, respectively (Table 6). Additionally, the determination coefficients between measured and predicted values varied from 65.0 to 96.0%. These results indicate the absence of constant and systematic errors, confirming the good reproducibility of the estimates from the proposed models when applied to the Caatinga vegetation.

Results of this study indicate that parametrization on daily and seasonal bases improved estimates for rain events where interception is less dependent on Sc than on Ec/R, that is, for larger and more intense events. During heavy storms the canopy tends to rapidly saturate, decreasing the influence of canopy storage capacity and increasing the control of Ec/R [8]. This behavior also explains the minor error propagation, observed when cumulative interception was simulated (Figure 5). However, for other ecosystems, this trend should be better investigated.

The sparse Gash model parametrized on a daily basis is indicated for vegetation densities that change gradually, but relatively rapidly, as well as for vegetation that changes more slowly, but is subject to infrequent rainfall [12]. It requires leaf area index monitoring and a more complex implementation, but represents important conceptual improvements in the rainfall interception simulations, giving accurate estimates from low to high intensity storms, as well as for events with different amounts. However, when it is not possible to use the more expensive instrumentation required for parametrizing this approach, or the greater data processing during simulations, the seasonal sparse Gash model is capable of considering the variability of Caatinga species regarding foliation and defoliation, which is reflected by the canopy and climate parameters associated with rainy and dry periods. It does not require leaf area monitoring and equations are simplified, resulting in a less complex simulation with reliable approximations.

#### *4.4. Limitations and Constraints*

The methodological challenges in measurement and data processing when modelling and validating rainfall interception are associated with the complex and expensive micrometeorological instrumentation, as well as the long data acquisition periods required in this process. Additionally, throughfall is highly spatially heterogeneous at small scales, while rainfall interception and stemflow are variable across species, requiring a measurement scheme capable of sufficiently take into account these differences, but also coherent with the financial and logistical constraints. For this, it is important to attempt to correctly locate the gauges and divide the area under each tree, also systematically performing the relocation of gauges by applying well defined methods and minimizing errors originating from spatial variability. These precautions were followed in this study, agreeing with other works [2,10,13,16,33] that had proved the possibility of monitoring throughfall, stemflow and rainfall interception by considering fewer sample trees per studied species in a credible way.

Among the Gash model parameters, canopy storage capacity and evaporation rate from wet canopy are the most difficult to obtain individually [21,47]. The quantification of canopy storage capacity can be determined experimentally for a particular species using laboratory methods [12]. However, the indirect methods are relatively low-cost and require no complex instrumentation, having been preferred in most of studies for Gash model parametrization [3,7,8,23], despite the long measurement period required. In this work, a regression-based method was applied, but future efforts for obtaining reliable measurements of this parameter or decomposing it into easily measurable physical components are encouraged [55]. The measurement of the evaporation rate from wet canopy also involves high costs and technical difficulties and this parameter is frequently estimated by means of the Penman–Monteith method [47,52]. This method was used in this

study, performing very well, despite requiring data of many micrometeorological variables. Other formulations and methods that overcome the high data input requirements of the Penman–Monteith equation have been proposed [2,24,34], but the best method for this purpose remains controversial and deserves further attention [47].

Finally, climate change is altering the water cycle and world ecosystems, with precipitation and phenological responses of plants to habitat being directly affected. Considering the Brazilian semi-arid region, the global climate change scenarios indicate that the aridity of the region will tend to increase in the next century, evidencing its vulnerability [18]. Therefore, a good comprehension of the hydrological processes and the study of their behaviors considering the climate change projections is essential for the evaluation of water availability and the anthropogenic effects in Caatinga. Thus, rainfall interception modelling based on these projections and the use of validated models for studies of land management are required.

#### **5. Conclusions**

Observed cumulative rainfall interceptions varied from 10.1 to 26.7% of gross rainfall during the rainy season and from 8.3 to 12.8% during the dry season for the five studied Caatinga species, indicating that significantly lower throughfall and stemflow reached the soil as available water input during the dry periods of 2016 and 2017. The frequency distributions of the event size and intensity did not differ significantly between dry and rainy seasons, with the very small (0.0–5.0 mm) and low-intensity (0–2.5 mm h−1) events being significantly more frequent. However, the highest percentages of total gross rainfall were verified for the very large events (30.1–40.0 mm) during the rainy season and for middle events (10.1–20.0 mm) when considering the dry season. The low intensity storms contributed to a higher percentage of gross rainfall only during the dry season, while middle and large intensity rainfall events (2.6–10 mm h<sup>−</sup>1) were responsible for the greatest percentages of gross rainfall during the rainy season. The sparse Gash model parametrized on a daily basis performed slightly better than the seasonal one, but both approaches resulted in very good or good agreement between the modelled and estimated interception, with cumulative mean relative errors between 1.23 and 8.98%. Seasonal model predictions were the most sensitive to canopy and climatic parameters, with canopy storage capacity presenting larger effects on the Caatinga rainfall partitioning, followed by Em, c and R, respectively, for both simulation approaches. Future works should focus on finding reliable methods for measuring canopy storage capacity, as well as on formulations capable of accurately estimating the evaporation rate from wet canopy requiring smaller input variables. Furthermore, the Gash model should be parametrized for other Caatinga species, with the validated approaches providing a basis for studies of land management, including the evaluation of degraded areas and effects of climate changes.

**Author Contributions:** Conceptualization, D.C.L. and A.J.S.N.; data curation, L.S.B.S. and C.A.A.S.; methodology, D.C.L., A.J.S.N., T.G.F.S. and S.Z.; resources, T.G.F.S. and S.Z.; validation, D.C.L. and A.J.S.N.; writing—original draft, D.C.L.; writing—review and editing, A.J.S.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Research Support Foundation of the Pernambuco State (FACEPE-APQ-0215-5.01/10 and FACEPE-APQ-1159-1.07/14), the National Council for Scientific and Technological Development (CNPq-475279/2010-7, 476372/2012-7, 305286/2015-3, 309421/2018- 7 and 152251/2018-9) and the Coordination for the Improvement of Higher Education Personnel (CAPES-Finance Code 001).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **Study of the Overflow Transport of the Nordic Sea**

**Wenqi Shi 1,\*, Ning Li 2,\* and Xianqing Lv 3,4**


**Abstract:** Changes in the climate system over recent decades have had profound impacts on the mean state and variability of ocean circulation, while the Nordic Sea overflow has remained stable in volume transport during the last two decades. The changes of the overflow flux depend on the pressure difference at the depth of the overflow outlet on both sides of the Greenland-Scotland Ridge (GSR). Combining satellite altimeter data and the reanalysis hydrological data, the analysis found that the barotropic pressure difference and baroclinic pressure difference on both sides of the GSR had a good negative correlation from 1993 to 2015. Both are caused by changes in the properties of the upper water, and the total pressure difference has no trend change. The weakening of deep convection can only change the temperature and salt structure of the Nordic Sea, but cannot reduce the mass of the water column. Therefore, the stable pressure difference drives a stable overflow. The overflow water storage in the Nordic Sea is decreasing, which may be caused by the reduction of local overflow water production and the constant overflow flux. When the upper interface of the overflow water body in the Nordic Sea is close to or below the outlet depth, the overflow is likely to greatly slow down or even experience a hiatus in the future, which deserves more attention.

**Keywords:** Nordic Sea; overflow flux; barotropic pressure; baroclinic pressure

#### **1. Introduction**

As an important driver of thermohaline circulation, the Nordic Sea overflow has a profound impact on environmental changes in the Arctic and even the world. In the Nordic Sea, high-density water bodies with a geopotential density (σө) greater than 27.8 kg/m<sup>3</sup> and shallower than the Greenland-Scotland Ridge (GSR) depth can overflow. There are three overflow channels on the GSR. From west to east, they are the Denmark Strait (DS), the Iceland-Faroe Ridge (IFR), and the Faroe-Shetland Channel (FSC) (Figure 1a,b). The overflow of dense water between Greenland and Shetland consists of the Faroe Bank Channel (FBC) overflow and Wyville Thomson Ridge (WTR) overflow, and FBC is the main channel for FSC overflow. The high-density water overflowing from these channels forms the North Atlantic Deep Water (NADW), which affects the nature of the deep-water mass and the deep circulation in the North Atlantic [1–4].

Theoretical analysis and field measurement results show that the Nordic Sea overflow is hydraulically controlled. In hydraulic control theory, changes of the overflow flux through a strait depend only on the total pressure difference at the depth of the sill on both sides of the GSR [1–4]. The total pressure difference is equal to the barotropic pressure difference plus the baroclinic pressure difference, depending on the local and remote physical processes, such as convection, mixing, and circulation, which further determine the overflow flux of the Nordic Sea [4].

**Citation:** Shi, W.; Li, N.; Lv, X. Study of the Overflow Transport of the Nordic Sea. *Water* **2021**, *13*, 2675. https://doi.org/10.3390/w13192675

Academic Editors: Andrzej Walega and Tamara Tokarczyk

Received: 21 July 2021 Accepted: 22 September 2021 Published: 27 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Figure 1.** (**a**) Topographic map of the Nordic Sea and (**b**) bottom depth map of the GSR. Note: The Nordic Sea is composed of the Norwegian Sea, the Icelandic Sea, and the Greenland Sea. In the Norwegian Sea, LB stands for the Lofoten Basin and NB stands for the Norwegian Basin; the DS, IFR, and FSC stand for the three main overflow channels in Greenland-Scotland Ridge, namely Denmark Strait, Iceland-Faroe Ridge, and Faroe-Shetland Channel, respectively; the purple square represents the selected point for calculating the pressure difference between the two sides of the GSR; the red line is the selected section for the bottom depth map. Bottom depth of the oceanic part of the GSR and overflow flux are based on Hansen and Østerhus [2].

Affected by climate change, the deep convection in the Nordic Sea has been weakened significantly from the 1960s to the beginning of the 21st century; by about 2006, the depth of deep convection in Greenland was less than 1000 m [5,6]. Recent studies showed that although there is a decreasing trend in atmospheric forcing from 1993 to 2016, the depth of convection in the Greenland Sea in winter has a deepening tendency. This is due to the increase in the salinity of seawater in the upper 1500 m, which results in the weakening of stratification inside the Greenland Sea circulation [7,8]. Modern climate models have found that the overflow flux of the Nordic Sea has a good consistency with the Greenland Sea convection, showing a weakening trend [9–12]. However, this weakening is not reflected by the measured data. The field measurement found that the overflow flux of the Nordic Sea remained strong and stable from 1995 to 2015, and there was no significant trend change [13,14].

For the phenomenon of stable overflow transport of the Nordic Sea during the last two decades, there have been studies explaining it from different aspects. Based on the model results, Olsen et al. [4] pointed out that the upper interface of the overflow water in the Nordic Sea declined from 1948 to 2005, which would cause a decrease in the pressure difference on both sides of the GSR. However, the rising sea level of the Nordic Sea offsets this effect, resulting in no trend change in the total pressure difference on both sides, making the overflow flux stable. Some other studies showed that the circulation of the Atlantic waters in the Nordic Sea has a greater impact on overflow changes, and the impact of weakened convection has been concealed [15]. Zhang and Thomas [16] believed that the Arctic Ocean, rather than the Greenland Sea, is the northern end of the mean Atlantic Meridional Overturning Circulation (AMOC). They further pointed out that the deep convection of the Labrador Sea and the Greenland Sea contribute the least to the mean AMOC, and AMOC may not be significantly weakened by the closure of the deep convections. However, some other studies still believed that the Greenland Sea is the main source area of the densest overflow water into the North Atlantic after 2005 and is the main ventilation area of the deepest layer in the North Atlantic [7,17]. Other studies suggest that the volume of the dense water above the GSR sill depth in the Nordic Seas is sufficient to supply decades of overflow transport without dense water production [1–3]. The premise in such estimations, however, is that all dense water above the sill depth is freely available for overflow transport. However, basin-scale oceanic circulation is nearly geostrophic and its streamlines are basically the same as the isobaths. The vast majority of the dense water is stored inside the closed geostrophic contours in the deep basin and thus is not freely available for overflow transport [18]. Therefore, an external force or a non-geostrophic mechanism is required to help transport the interior water mass to the boundary current. The numerical simulation results of Yang and Pratt [19] show that 80%–85% of the dense water above the GSR sill depth in the Nordic Seas is not freely available for overflow transport, and the amount of the dense water freely available to overflow accounts for only 15%–20%. Therefore, the Nordic Seas has a relatively small capacity as a dense water reservoir and thus the overflow transport is sensitive to climate changes.

In short, there is still controversy about the reasons for the stable overflow flux in the past two decades. Based on the satellite altimeter data and the reanalysis hydrological data, this paper will analyze the changes in the barotropic pressure and baroclinic pressure on both sides of the GSR and then discuss the reasons for the long-term stable flux of the Nordic Sea overflow by the hydraulic control theory.

The structure of this paper is as follows. Chapter 2 introduces the data and the method for calculating the pressure. Chapter 3 evaluates the credibility of the EN4 data to calculate the pressure by comparing the measured hydrological data and the overflow flux results. Chapter 4 mainly analyzes the spatial distribution of the change trends of the positive pressure, baroclinic pressure, and total pressure on both sides; the change characteristics of the pressure difference on both sides of the GSR; the changes in depth of the overflow water interface in the Nordic Sea; and then analyzes the reasons for stable overflow flux from 1993 to 2015. Chapter 5 mainly analyzes the correlation between the positive pressure and baroclinic pressure on both sides of the GSR and the role of the changes in the properties of the upper seawater. Chapter 6 is the conclusion.

#### **2. Data and Methodology**

#### *2.1. Satellite Altimeter Data*

The Sea Level Anomaly (SLA) data in this paper is monthly averaged data from 1993 to present of the DUACS 2014 database from the French Space Research Center (CNES), which is merged with multi-satellite altimetry data (Available online: http://www.aviso. oceanobs.com/duacs/ (accessed on 5 May 2016)). The data use Mercator projection with a horizontal resolution of 1/4◦ × 1/4◦, and are corrected by atmospheric pressure correction,

tide correction, and dry tropospheric correction. The SLA data of DUACS 2014 are based on the Mean Sea Surface (MSS) from 1993 to 2012. Since the mean sea level is the height on a fixed earth reference ellipsoid, the SLA contains sea level change signals caused by the relative crustal movement during this period, mainly as the Glacial Isostatic Adjustment (GIA). Tamisiea and Mitrovica [20] gave the distribution map of the GIA effect on the sea level change measured by the altimeter (their Figure 3b), and their results showed that the GIA has an effect of no more than 0.15 mm/yr on the sea level change trend in the Nordic Sea, so it can be ignored here. Volkov and Pujol [21] verified through field measurement that AVISO's altimeter data can be used to study sea level changes and surface currents in the Nordic Sea.

#### *2.2. Hydrological Data*

The hydrological data in this paper are the monthly averaged reanalysis data of the EN4 hydrological data set of the Met Office [22]. The data are obtained from a large amount of observational data, mainly including WOD data (World Ocean Database), GTSPP data (Global Temperature and Salinity Profile Project), Argo data, and ASBO (Arctic Synoptic Basin-wide Observations) data, which have high credibility. The horizontal resolution of the data is 1◦ × 1◦, and the coverage area is 1◦ E–360◦ E and 83◦ S–89◦ N. The data has 42 vertical layers, and the thickness of the water layer ranges from 10 m in the upper layer to 300 m in the deep ocean.

Since the grid of SLA is inconsistent with the hydrological data grid, the SLA data needs to be interpolated to the same grid point as the hydrological data. The interpolation method is used to obtain the mean value of the 16 SLA grid points in each hydrological grid. If there are more than 8 missing values in a SLA data grid, the mean value at that point is also assigned as missing.

The hydrological observation data of the "Mike" station (Ocean Weather Ship Station Mike, here referred to as OWS-M station) comes from the European Ocean Observatory Network (Euro SITES, Available online: http://www.eurosites.info/stationm.php (accessed on 26 December 2016)). The station is located in the center of the Norwegian Sea (66◦ N, 02◦ E) and provided long-term ocean and meteorological profile data almost daily from October 1948 to November 2009.

#### *2.3. Method for Calculating Pressure*

Based on the hydrostatic assumption, the pressure at a certain depth without considering the atmospheric pressure is:

$$\mathbf{P} = \mathbf{P}\_{\text{trop}} + \mathbf{P}\_{\text{clin}} = \rho\_0 \mathbf{g} \boldsymbol{\zeta} + \mathbf{g} \int\_{\mathbf{z}}^{0} \mathbf{p} \, \mathbf{z} \tag{1}$$

where Ptrop represents the barotropic pressure and Pclin is the baroclinic pressure; ρ<sup>0</sup> = 1028 kg/m<sup>3</sup> is the surface seawater density; g = 9.8 m/s<sup>2</sup> is the gravitational acceleration; ζ is the sea level height, and here is taken as SLA; ρ is the seawater density, which is derived from the temperature and salt data of EN4; z is the calculated pressure depth, and unless otherwise specified, it is taken as 840 m, which is the maximum depth of the GSR sill. Actually, the mean sea surface level is not horizontal, and the spatial difference is huge. However, since this article focuses on the temporal change of pressure rather than the absolute value, taking ζ as the Sea Level Anomaly (SLA) will not affect the final analysis result.

#### **3. Applicability Analysis of EN4 Data**

#### *3.1. Comparison with Observations at OWS-M Station*

Analysis of the data from the OWS-M station shows that there are enough data for Pclin calculations every month above 1000 m depth. The monthly and annual mean results of Pclin calculated from the data are shown in Figure 2. Comparing the results of EN4 data with the results of EN4 data, their annual mean change curves have a high degree of overlap. The decline rate of annual mean Pclin at the OWS-M station during 1949–2009 was −0.55 ± 0.26 × <sup>10</sup><sup>2</sup> Pa/dec, while the EN4 data was −0.3 ± 0.26 × 102 Pa/dec, consistently showing a downward trend. Considering the annual mean Pclin at the OWS-M station had some missing data, the difference in the decline rate between the two data sets may have been caused by the missing data.

**Figure 2.** The Pclin anomaly from OWS-M data and EN4 data. Note: The OWS-M station data are the daily profile observation results. Firstly, the monthly mean temperature and salinity values at different depths are obtained by averaging, and then the monthly mean Pclin (at 840 m depth) is calculated from temperature and salinity; when the cumulative observation level of a month is less than 10 layers or the observation depth is less than 80 data 0 m, the month will be treated as a missing measurement. The annual mean Pclin of the OWS-M station is obtained from the monthly data average. A year is regarded as a missing year if there are more than 4 months in the year of missing annual mean data and is not shown.

#### *3.2. Compared with the Observed Overflow Transport*

Based on mooring ADCP and temperature and salinity observations, Hansen and Østerhus [3] found that there was no significant trend change in FBC overflow flux from 1995 to 2005, and the trend change did not exceed 0.2 Sv, which is only 10% of the mean flow. Figure 3a–c shows that the SLAs on both sides of the GSR are both increased during this period, while the Pclin at the depth of 840 m is decreased at the same time. The final P obtained has no remarkable trend change in the Norwegian Sea or in the Icelandic Sea. Therefore, the trend changes characteristics of the FBC overflow obtained from SLA and EN4 data are more credible.

The measured data show that the trend of DS overflow flux in the 15-year period from 1996 to 2011 is −0.4 Sv. However, the trend is below the 70% confidence level of the *t*-Test, so it is not significant [23]. Here the pressure at the depth of 640 m (approximately the depth of DS) on the north and south sides of the DS increases by the same magnitude, and the pressure difference between the two sides is basically unchanged. The spatial distributions of Pclin and P at the depth of 640 m are basically the same as in Figure 2a–c, which is not shown separately here. Thus, the calculation results here are consistent with the observation.

**Figure 3.** The change trend of SLA (**a**), Pclin (**b**), and P (**c**) from 1995 to 2005. The gray area represents the sea area that has not passed the 95% significance test; the white area represents the sea area where the period of missing data is longer than 7 years; the black lines are the 1000 and 3000 m isobaths. The pressure is calculated based on a reference depth of 840 m.

#### **4. Change from 1993 to 2015**

In the published literatures, the observation of GSR overflow flux is available until 2015 [14]. Since the Nordic Sea overflow is hydraulically controlled, the pressure difference on both sides of the GSR can be used to analyze the long-term changes of the overflow flux. The depth of the deepest GSR sill is about 840 m on the FBC, which can be used to calculate the pressure difference [4]. As shown in Figure 4, from 1993 to 2015 the SLAs of the Nordic Sea and the North Atlantic subpolar region near the GSR basically increased at the same rate; Pclin mostly declined in the south of GSR, increased near DS in the north of GSR, and had no significant change in the south of Norwegian Basin. The pressure difference in the west of Iceland had a clear upward trend; the pressure difference to the east of Iceland was basically unchanged. This means that DS overflow increased, while FBC overflow and IFR overflow did not change much. Therefore, the total overflow in the Nordic Sea slightly increased.

**Figure 4.** SLA (**a**), Pclin\_840 (**b**), P\_840 (**c**), P\_bottom (**d**) changes trends from 1993 to 2015. The gray area represents the sea area that has not passed the 95% significance test; the white area represents the sea area where the period of missing data is longer than 7 years. The black lines are the 1000 and 3000 m isobaths; the purple square represents the selected point for calculating the pressure difference between the two sides of the GSR in panel (**c**).

In the Nordic Sea, the depth of σ<sup>ө</sup> = 27.8 kg/m3 (D27.8) (Figure 5a,b, the upper interface of the overflow water) is more consistent with the spatial distribution of the change rate of Pclin, indicating that changes in the properties of the upper seawater necessarily indicate the adjustments of the upper interface of the overflow water. Especially in the Nordic Sea, the sinking of D27.8 in the Lofoten Basin is about 100 m/dec, which may be directly caused by the reduction of deep convection in the Greenland Sea [5,6] or the weakening of other dense water production. When the total overflow transport flux remains unchanged, the amount of overflow water flowing out of the Lofoten Basin almost remains unchanged. Therefore, the reduction of dense water supply leads to the rapid sinking of the overflow water interface in the Lofoten Basin. There, the rapid decline of Pclin and the rapid rise of Ptrop occur at the same time, while P is basically unchanged, indicating that the above changes are probably caused by the change of physical property of the upper water.

**Figure 5.** The spatial distribution of the mean value (**a**) and change rate (**b**) of D27.8 from 1993 to 2015. The D27.8 is the depth of σө= 27.8 kg/m3.

The Labrador Sea is a fast-sinking area of D27.8, with a sinking rate up to 160 m/dec or more, which is consistent with the reported weakening of convection there [24,25]. In the modern climate, the Nordic Sea overflow, entrainment process, and Labrador Sea convection provides about 1/3 of the deep branch of the radial overturning circulation in the Atlantic Ocean [4,26]. Under the condition that the overflow of the Nordic Sea is relatively stable, the convection in the Labrador Sea is significantly weakened, which may be the main reason for the significant weakening of the AMOC near 25◦ N [24,27]. However, some relatively new observational evidence has indicated that the deep convection of the Labrador Sea has the smallest total contribution to the subpolar overturning circulation [28,29].

The convection of the Labrador Sea is significantly weakened, which causes the upper interface of the dense water to sink quickly. Since the upper interface of the overflow water in the Labrador Sea is deeper than 1500 m, its impact on the Pclin at 840 m is small, and the decrease rate of Pclin at 840 m depth is only −<sup>2</sup> × 102 Pa/dec. The deepening of the overflow water in the Labrador Sea means that warming and freshening of the entire water column causes a large increase in SLA, which is consistent with the calculated results. The greater the bottom depth is, the greater the increase of SLA is. However, there is no significant trend change in the mass of the entire water column from surface to the bottom (Figure 4d).

Although the changes in properties of seawater can ensure the mass conservation of the whole water column, the compression or expansion of the water column caused by the change of properties of seawater will lead to the change of the mass ratio of the upper and lower water column at a certain depth. Therefore, in hydrostatic balance, the pressure change of seawater at a certain depth may be caused by the change of properties of seawater below this depth, and the change of properties of seawater above this depth has no effect on it. The different changing trends of the pressure at 840 m depth and the seabed depth in the Labrador Sea and the Irminge Sea in the south of Greenland (Figure 4b,d) show this effect.

Two points have been selected at the upstream and downstream ends of the FBC to construct the temporal variations of pressure difference. Based on the overflow water sources in different overflow channels and combining the location given by Olsen et al. [4], we selected (64 N, 2 W) and (62 N, 15 W) to estimate the FBC overflow flux (the location is shown in Figure 1). It can be seen from Figure 6 that the inter-annual variation characteristics of ΔPtrop and ΔPclin obtained in this paper are very consistent with Figure 2 of Olsen et al. [4]. Both results showed the minimum values of ΔPclin and ΔP in 1995, and

the relative maximum values of ΔPtrop and ΔP in 2003; from 1993 to 2005, ΔPclin and ΔP were rapidly rising and ΔPtrop had no significant change trend. At the same time, the inter-annual variation of <sup>Δ</sup>P calculated here is about 2 × 102 Pa and about 10% of the mean ΔP, which is basically consistent with the observed inter-annual variation of FBC overflow [3]. In short, the pressure difference between the two points selected in this paper can be used to estimate the FBC overflow flux change. From the spatial distribution map of the SLA trend rate (Figure 4), it can be seen that the trend rate of SLA has a good spatial continuity in the sea areas near both sides of the GSR, and the results would not be significantly changed due to slight difference in the selection of the grid location.

**Figure 6.** The variations of annual mean of ΔPtrop, ΔPclin, and ΔP. The positions of the selected grid points to calculate pressure difference on both sides of the GSR are shown in Figure 1.

Specifically, ΔPtrop experienced a slow decline with fluctuation from 1993 to 2005 and reached the minimum value in the past 23 years before 2005. ΔPtrop increased with fluctuation from 2005 to 2013 and rose rapidly from 2013 to 2014; after that, it fell back. The year 2014 had the maximum value of ΔPtrop in the past 23 years (Figure 6). ΔPclin first decreased slightly in the period of 1993–1997, then increased before 2004, and reached the maximum value in the past 23 years before 2004; then it decreased slowly in the period of 2004–2013, and decreased rapidly in the last two years. ΔPclin in 2015 reached the minimum value in the past 23 years. ΔP was basically at an average level in 1993, followed by a relatively large fluctuation. After experiencing the minimum value in 1995 and the maximum value in 2003, it basically returned to the mean level in 2015. The linear regression of the annual mean <sup>Δ</sup>P results in a change rate of 1.6 × <sup>10</sup><sup>2</sup> Pa/dec. Olsen et al. [4] gave a linear coefficient of FBC overflow flux change (Δq) and ΔP of k = 10−<sup>3</sup> Sv/Pa. Using this linear coefficient, we obtain the FBC overflow enhancement rate of about 0.16 Sv/dec, which is quite small relative to the mean FBC overflow flux (2.9 Sv). At the same time, the linearly increasing trend of ΔP failed the 95% confidence test but passed the 90% confidence test. In fact, <sup>Δ</sup>P in 2015 was only about 1 × <sup>10</sup><sup>2</sup> Pa larger than in 1993, which is quite small.

The changes in these three parameters have no significant correspondence with NAO, and most of the wind stress curl changes in the Nordic Sea are related to NAO [30]. This indicates that the interannual sea level changes are not mainly driven by wind stress, but more likely are the result of changes in the properties of the upper seawater.

#### **5. Relationship between Barotropic Pressure and Baroclinic Pressure**

Olsen et al. [4] concluded that ΔPtrop and ΔPclin on both sides of the FBC have a correlation lag of about three years, and analyzed the mechanism of the correlation as follows: due to wind stress, the sea level difference on both sides of the FBC increases (decreases) and the ΔP on both sides increases (decreases) through the barotropic pressure effect. Then, the overflow transport is enhanced (weakened), causing the iso-density interface in the Norwegian Basin to sink (rise) and the ΔPclin gradually decreases (increases); and then ΔP gradually decreases (increases) until recovers to normal level. This feedback mechanism could help ΔP remain stable, which means the overflow transport is stable. They use a simplified two-layer model to express the mechanism as:

$$\mathbf{P}\_{\text{trop}} = \rho\_0 \mathbf{g} \Delta \mathbf{h} \tag{2}$$

$$
\Delta \mathbf{P}\_{\rm clin} = -\mathbf{g} \Delta \boldsymbol{\varrho} \Delta \mathbf{D} \tag{3}
$$

$$
\Delta \mathbf{P}\_{\rm clin} = -\mathbf{g} \Delta \rho \Delta \mathbf{D} = \frac{-\mathbf{g} \Delta \rho \mathbf{k} \mathbf{T}}{\mathbf{A}} \Delta \mathbf{P}\_{\rm top} \tag{4}
$$

where <sup>ρ</sup><sup>0</sup> = 1.025 × 103 kg/m3 is the surface seawater density, g = 9.8 kg/m3 is the gravitational acceleration, and Δρ = 0.5 kg/m<sup>3</sup> is the density difference between overflow water and upper water body. Linear regression coefficient of overflow flux change (Δq) and pressure difference (ΔP) is k = 10−<sup>3</sup> Sv/Pa. A is the contact area between the overflow layer and the upper layer in the Nordic Sea, or rather the area of the Norwegian Sea deeper than 500 m, which is equal to 5.8 × 1011 <sup>m</sup><sup>2</sup> [4]. T is the time for the high-density water interface to sink ΔD after the barotropic pressure disturbance, which is also the time for ΔP to restore to the initial state. The calculated T is approximately equal to three years.

The monthly mean variation of ΔPtrop, ΔPclin, and ΔP was constructed based on EN4 and SLA data, and the correlation between ΔPtrop and ΔPclin lagging or leading in different months was analyzed (Figure 7). When ΔPtrop is about three months ahead of ΔPclin, the negative correlation between them is the largest (−0.59). Olsen et al. [4] defined the horizontal spatial area occupied by overflow water as the seabed deeper than 500 m. However, the dense water in the Norwegian Sea is not freely available for overflow transport, and the dense water in the center of the basin, which occupies most of the area, is circulated by the boundary oceanic circulation. Therefore, the size of the effective overflow area is much smaller than that of the ocean basin. Based on the feedback mechanism of Olsen et al. [4] and the lag time obtained in this article, the horizontal spatial range of available overflow water upstream of the FBC can be estimated to be 0.5 × 1011 <sup>m</sup>2, which is about 1/12 of the value given by Olsen et al. [4]. This ratio is close to the percentage of available overflow water in the total overflow water in the Nordic Sea (80%~85%) obtained by other studies [19].

**Figure 7.** The variation (**a**) and the correlation for different lag time lengths (**b**) of ΔPtrop and ΔPclin.

The high Pclin correlation with the station in the southern part of the Norwegian Sea (64 N, 2 W) is limited to a small area in the southern part of the Norwegian Sea (the area with a correlation greater than 0.8 in Figure 8). The Pclin in this area has a good consistency of change, which can be considered as the available overflow area upstream of the FBC overflow. The area with a correlation greater than 0.8 is about 1.8 × <sup>10</sup><sup>11</sup> m2, and the area with a correlation greater than 0.9 is 0.9 × <sup>10</sup><sup>11</sup> m2. It is more likely that the changes of ΔPtrop and ΔPclin are both dominated by changes in seawater properties, so the largest negative correlation between them basically has no lead or lag (Figure 9).

**Figure 8.** The correlation between Pclin and Pclin at selected stations in the Norwegian Sea. The time series of the correlation analysis has undergone a 12-month moving average processing. The correlation here is the maximum correlation within 5 years of lead or lag time. The selected stations in the Norwegian Sea are shown in Figure 1 as a purple square. The time series of ΔPtrop, ΔPclin, and ΔP are all carried on a 12-month moving mean to remove seasonal fluctuations in this figure. For the sake of comparison, ΔPclin and ΔP are shown as anomalies.

**Figure 9.** The correlation between Pclin and Ptrop.

Changes in the properties of the upper seawater will cause the reverse change of Pclin and Ptrop, while the total pressure will not change due to the unchanged seawater quality. Therefore, the sea area with a stronger negative correlation between Pclin and Ptrop indicates that changes in the properties of the upper seawater play a greater role in the changes of both there. It can be seen from Figure 9 that there is a strong negative correlation between Pclin and Ptrop in the southern sea area of GSR, the correlation coefficient is close to −1.0, while the total pressure at this place has no trend change characteristics (Figure 4), which shows that the changes of Pclin and Ptrop are mainly caused by the changes in the properties of the upper seawater. In the Nordic Sea north of GSR, this negative correlation is not so strong. Among them, in the Norwegian Sea, Pclin and Ptrop have a certain negative correlation, indicating that the change in the properties of the upper seawater is one of the important factors which cause the changes in the two. There are other processes that lead to the increase in the quality of the upper seawater, which causes a slight increase trend in the total pressure (Figure 4). The negative correlation between Pclin and Ptrop is no longer significant in other areas of the Nordic Sea except the Norwegian Sea. In the Icelandic Sea, the correlation between Pclin and Ptrop is poor and the SLA increases significantly, which leads to a significant increase in the total pressure (Figure 4). There is a weak positive correlation between the two in the Greenland Sea, indicating that the changes of Pclin and Ptrop in this area are mainly affected by other processes.

Under hydrostatic assumption, changes in the density of seawater above 840 m depth will not change the hydrostatic pressure at this depth. To change the pressure at this depth, it needs to change the mass of the water column at this depth. There are two ways. One is to change the absolute mass of the water column, or to change the sea level through wind stress curl, runoff input, sea-air material flux, and other factors. The other is to change the relative mass of the water column by changing the density of the deep layer, causing the column to expand or contract. The mass percentage of the water column above the 840 m depth can change the entire water column.

At present, most ocean numerical models are based on Boussinesq approximation, which cannot reflect sea level changes caused by changes in seawater properties. When the density of the sea layer in Northern Europe decreases, the pressure obtained by simulation decreases, which in turn leads to the weakening of simulated overflow [31]. It can be seen from the results of this paper that the steric effect contributes to most of the sea level trend changes in the sea area surrounding the GSR and has a significant impact on the long-term changes in overflow transport. Therefore, the simulation and prediction of long-term changes in overflow requires the use of non-Boussinesq ocean models, considering the impact of changes in seawater properties on SLA.

#### **6. Conclusions**

The Nordic Sea overflow is hydraulically controlled; the changes of the overflow flux depend only on the pressure difference at the depth of the overflow outlet on both sides of the GSR. Based on the satellite altimeter data and the reanalysis hydrological data, we obtained a slight increase in the pressure difference between the two sides of the GSR from 1995 to 2015. However, this trend is not significant and is more consistent with the observed stable overflow flux. Among them, the barotropic pressure and baroclinic pressure in the southern sea area of the GSR have a very good negative correlation (correlation coefficient is close to −1.0). The changes in both are basically caused by the changes in the properties of the upper seawater, and the total pressure there is only a slight increasing trend. The barotropic pressure and baroclinic pressure of the Norwegian Sea in the northern part of the GSR have a certain negative correlation (correlation coefficient is about −0.6), indicating that changes in the properties of the upper seawater are important factors that cause changes in the barotropic and baroclinic pressures in the sea area, and other processes can also lead to a slight increase in the barotropic pressure there. While the correlation between the barotropic pressure and the barotropic pressure in the Icelandic Sea is poor, the barotropic pressure increases significantly which leads to a significant increase in the total pressure there.

By selecting two representative points, the barotropic pressure difference and baroclinic pressure difference on both sides of the FBC are constructed. The changes in the barotropic pressure and baroclinic pressure on both sides of the FBC are more likely caused by the changes in the properties of the local upper seawater. The total pressure difference caused no significant trend changes characteristics between 1993–2015, which is consistent with the observation of stable overflow flux.

In the Nordic Sea, the area with the fastest sinking of the overflow water upper interface is the Lofton Basin, with a sinking speed of more than 100 m/dec, indicating that the storage of overflow water there is rapidly decreasing. The physical processes that produce dense water, such as deep convection in the Greenland Sea, are weakening, and the source of overflow provided is reducing, leading to warming and lightening of the upper layer of the Norwegian Sea and sinking of the upper interface of the overflow water. However, the changes in the properties of the upper seawater in the Norwegian Sea cannot reduce upstream pressure in the depth of the sill to weaken overflow transport. Therefore, it will cause the upper interface of upstream overflow water to further decrease. In the future, when the depth of the overflow water upper interface in the Nordic Sea is less than the depth of the sill on the GSR, the overflow may greatly slow down or even experience a hiatus. This is worthy of close attention and further study.

**Author Contributions:** Conceptualization, W.S.; methodology, W.S. and N.L.; software, W.S. and X.L.; formal analysis, W.S.; writing—original draft preparation, W.S.; writing—review and editing, N.L.; supervision, X.L.; Funding acquisition, W.S. and X.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant number 41806219 and 41806002, the National key research and development program, grant number 2018YFC1407602 and 2019YFC1408405, and the Youth Science and Technology Star Project of Dalian (2020RQ020).

**Institutional Review Board Statement:** The study did not involve humans or animals.

**Informed Consent Statement:** We choose to exclude this statement as the study did not involve humans.

**Data Availability Statement:** The sea level anomaly (SLA) data in this paper are monthly averaged data from 1993 to present of the DUACS 2014 database from the French Space Research Center (CNES), which is merged with multi-satellite altimetry data (http://www.aviso.oceanobs.com/ duacs/ accessed on 24 September 2021). The monthly mean reanalysis data of the EN4 hydrological data set are from Met Office Hadley Centre observations data sets (https://www.metoffice.gov.uk/ hadobs/en4/ accessed on 24 September 2021).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

