Characterizing estuarine salinity patterns with event duration and frequency of reoccurrence approaches

Salinity is widely recognized as a dominant environmental variable in estuarine systems. However, drawing clear causal relationships between salinity and estuarine floral or faunal patterns of abundance, population status, or vigor is often confounded by the great spatial and temporal variation exhibited by salinity. Averages, or other statistics, over monthly, seasonal, or longer periods are often used, but may fail to capture patterns of ecological significance. This study introduces a novel approach to characterize estuarine salinity based on deriving events of specific magnitude and duration from a long‐term record; events that are also amenable to analysis of the frequency that they reoccur. With sufficient point data for salinity over the spatial domain of an estuary, these variables can be mapped and contoured to reveal spatial patterns of potential ecological significance. Model‐predicted salinity for a large estuary system of the Texas coast was used to develop and test these new techniques. While this system has great variation in salinity, both temporally and spatially, we find that this variability can be characterized into well‐defined and reoccurring salinity‐based events. Furthermore, the frequency of reoccurrence of those events over a long‐term record reveals patterns not captured by more common techniques such as averaging over similar time frames. These novel techniques provide an integrated spatial‐temporal approach to salinity pattern portrayal which may be useful for examining short‐term (daily–monthly) estuary processes and longer‐term (yearly–decadal) limits to estuary species distributions and habitat composition.

Estuaries, the transition zones between freshwater and saline marine environments, are fundamentally characterized by spatial variation in salinity. Considerable temporal variation is also exhibited due to changes in freshwater inflow volumes, tidal forces, wind, and evaporation, among other drivers (Ward and Montague 1996). This amalgam of spatial and temporal salinity conditions coupled with the variability in tolerance of individual species and communities, leads to the recognition of salinity as perhaps the most significant variable in an estuarine ecosystem (Alber 2002). In spite of its importance, methods of characterizing this variable are currently limited and efforts to better integrate and portray this spatial and temporal variation may lead to more thorough understanding of estuarine processes and conditions. The roles that salinity and its variation can play in an estuary are nearly countless in number, but one potential classification of these would utilize two realms: a habitat/structure focus and a process/function emphasis. The habitat/structure focus examines spatial patterns of salinity that reoccur with some regularity and are viewed as a primary determinant of macro-habitat types, location, composition, and condition (e.g., Hopkins et al. 1973;Howard and Mendelssohn 1999;Doering et al. 2002;Weilhoefer et al. 2013;Borgnis and Boyer 2016). The process/functional realm focuses primarily on salinity and its variation as a driver of various biological and ecological processes such as growth, production, spawning, disease progression, predation, and migration (e.g., Cain 1973;Drinkwater and Frank 1994;Thorpe 1994;Soniat et al. 2012). There are also salinity-modulated physical processes, such as density-induced circulation, independent of tides, and other forces (Pritchard 1952).
In these manifold investigations, a wide variety of temporal scales are utilized. Generally, studies of salinity's influence on fixed habitats and sessile community composition tend to focus on multiyear or seasonal time scales (e.g., Wetzel and Kitchens 2007;Borgnis and Boyer 2016). Seasonal time scales are also usually embedded in the habitat suitability paradigm (HSP) in widespread use to inform freshwater management approaches based on salinity preferences or needs for individual species (e.g., Cake 1983;Clark et al. 2004;Stachelek and Dunton 2013). Studies of "disturbance" with intense salinity changes over short time scales driving quick habitat instability or intense upsets to normal processes, tend to focus on exceptional events on the order of weeks or months (e.g., Jiang et al. 2014;Kinney et al. 2014). Research on the role of salinity and salinity dynamics in modulating estuarine biological and ecological processes (spawning, growth, seed production, disease progression, migration) tend to utilize time scales ranging from seasonal (e.g., Powell et al. 2003) down to daily (e.g., Cain 1973;Stazisar et al. 2015).
Despite the strong influence that salinity and/or salinity variation exert on estuarine habitats and biota, and the various time scales of importance, there are limited methods currently available for portrayal and analyses of this environmental parameter. Due to inherent spatial and temporal variation of salinity, portrayals typically fall into just two distinct categories: (1) a spatial (map) view, or (2) a time series (temporal) depiction at a limited number of physical locations.
The spatial approach, by necessity, generally adopts either a point-in-time ("snapshot") view, or, at the other end of the spectrum, a time-averaged view over a defined period such as specific days, months, seasons, or years. This is the common map approach as illustrated in Fig. 1 with a long-term average salinity portrayal. While a map's inherent focus on spatial variation is an essential tool, this time-averaged approach can prove limiting in explanatory power since it does not readily lend itself to tracking any underlying timedynamic behavior of salinity. In contrast, the time-explicit approach, usually based on a series of salinities, focuses on variation through time, but typically at a single or very limited number of points. This approach obviously brings salinity dynamism (or stasis) into sharp focus, but at the expense of a reduction in spatial context.
Herein, we propose a novel method that integrates both spatial and time-based approaches. For display purposes, we utilize a spatial portrayal, but we overcome the typically static view of that approach by depicting variables that explicitly embody temporal salinity behavior. The timebased variable depicted on the map can range from short periods with salinity in an exceptional range (e.g., 30 consecutive days with salinity above 20 in a typically oligohaline [0.5-5] zone) to long periods of relative stasis (e.g., 120 consecutive days of salinity in the 10-15 range). Additionally, our method can spatially depict the frequency with which these "events" may reoccur over long time frames, providing a framework for understanding whether any single event is Fig. 1. The Guadalupe and Mission-Aransas Estuary system located along the central Texas coast and the principal river sources of freshwater. Also shown are the average salinity values throughout the system for the period 1987-2009. The specific reference nodes A through G will be used throughout.

Johns and Heger
Salinity duration and frequency patterns common (e.g., once per 1 yr) or rare (e.g., once per 20 yr). We believe these portrayal methods will provide additional insights on how salinity may influence estuarine habitats and ecological processes.

Materials and procedures
This new approach develops two new variables derived from underlying salinity data, either measured or synthetic (modeled or otherwise predicted), to describe salinity patterns. We first introduce the variable consecutive salinity days (CSD) to represent the consecutive days in which salinity is within a certain fixed range at a given point within the estuary. Since this variable has a time scale of days, the underlying salinity data for derivation should be similar: at the daily average, or tidal cycle average or other single values per day (e.g., fixed time of day each day). The CSD variable integrates salinity magnitude (e.g., ranging from 10 to 20) and duration of consecutive occurrence (e.g., 120 d or longer) as a salinity "event" at a particular spatial location. The computed CSD reveals point-specific information regarding salinity behavior over a given time frame, such as the maximum duration within a specific year that salinity was in a specified range. These data can then be mapped if enough point data for CSD are developed (examples in the "Assessment" section).
The focus on consecutive days is utilized due to the cumulative effect that salinity stress, from either elevated or depressed values, often exerts on biota. For example, infection of Eastern oysters by Perkinsus marinus, is a cumulative pathogen cell division process when salinity is elevated above 20 and temperature above 208C (Hofmann et al. 1995). A cumulative period focus is also consistent with other environmental stressors such as periods of low dissolved oxygen (e.g., Best et al. 2007). Other nonconsecutive formulations such as that of Huang (2010) integrate salinity magnitude and temporal frequency of occurrence, but are not amenable to categorizing salinity behavior as continuous "events." The second new pattern-depiction variable, return period of consecutive salinity days (rCSD), measured in years, characterizes the maximums of the yearly CSD events by their frequency of reoccurrence (e.g., reoccurring at least once per 2 yr) over a long-term record. Thus, the reoccurrence variable, rCSD, can characterize any given salinity "event" as to how rare or common it is, based on its frequency of occurrence in the historical record. We will follow the format of other disciplines such as hydrology (Chow et al. 1988) or meteorology (Hershfield 1961) in portraying the frequency data by referring to the return period, the expected value of the interval between reoccurrence of the same events (e.g., 5 yr). There is an inverse relationship between return periods and frequency of reoccurrence. For instance, very common events would be those that typically reoccur with a return period interval of 1 yr or 2 yr while rare events would be typified by a return period of 20 yr.
While the underlying salinity data for deriving these new variables may be actual measured values, a typical monitoring network will likely be too sparse or have too short of a period of record to be suitable for this method's approach of developing a map. Thus, the measured data will almost certainly need to be supplemented with data synthesized with statistical or simulation methods as utilized herein. These techniques may also be required to fill in missing temporal data at actual measurement sites.
To illustrate this method, we will rely on a set of modelsimulated salinity data for the interlinked Guadalupe and Mission-Aransas Estuary system, which covers approximately 1025 km 2 on the central Texas coast (Fig. 1). This estuary system exhibits pronounced variability in freshwater inflows and salinity, both intra-annually and over multiyear time scales, the later driven largely by El Niño Southern Oscillation periodicity (Tolan 2007). The great spatial and temporal variation in salinity in this estuary system provides an excellent setting to illustrate our approach. Figure 2 illustrates the range of salinity at several locations in the estuary system over the period 1987-2009.
The method developed and presented here is greatly aided by the fact that the Texas Water Development Board (TWDB) maintains a calibrated mathematical model, known as TxBLEND, which predicts the vertically averaged (twodimensional) salinity within this shallow estuary system based on freshwater inflows and other independent variables (e.g., tides and winds). This model was recently recalibrated

Johns and Heger
Salinity duration and frequency patterns ) and updated to include revised inflow estimates . The period of record covered by TxBLEND for this study was 01 January 1987-31 December 2009, or 23 yr. TxBLEND subdivides the estuary into a fine mesh of interconnected nodes (approximately 1900 nodes in the contoured area of Fig. 1) and simulates the salinity at each with a time step of 3 min with output generated at 1 h intervals. Important for this study is that TxBLEND provides a fine spatial scale which facilitates the derivation of this study's salinity variables. We utilized model-generated daily average values of salinity at a set of 162 well-dispersed model nodes (Fig. 3). This level of resolution was found to be adequate to cover the entirety of the interlinked Guadalupe and Mission-Aransas Estuary system and provide suitable spatial density for the contouring and mapping exercises for CSD and rCSD discussed below. Of those selected nodes, 147 were employed for spatial interpolation and contour map development and 15 were reserved for validation purposes discussed later. A subset of nodes, called reference nodes on Fig. 3, will be used throughout to illustrate important methodological concepts and results. Figure 4 illustrates an example of the TxBLEND model's time-series prediction of daily salinity for the 2004-2005 period at two reference nodes within the Guadalupe Estuary (as presented in Fig. 1). At both locations, there is a clear response of lowered salinity during periods of high to very high inflows and the opposite response during low inflow periods. For example, in May through June 2004, average daily inflows were approximately 23 million m 3 d 21 and were sufficient to lower salinity at reference node A to 0-1 for much of this period while at B salinity ranged from 0 to 5. Another even higher inflow event, with an average inflow of approximately 66 million m 3 d 21 for the month of November forced salinity at both locations to zero for most of the month. Except for these extreme inflow magnitude events, however, at reference node B, salinity is consistently at a higher level due to its location farther down the estuary away from the Guadalupe River, the major freshwater source.

Salinity duration events
As introduced above, the variable CSD 0-10 denotes a count of consecutive days in which salinity is within the range of 0-10 at some specific fixed point in the estuary system. Figure 5 illustrates the process for determination of this variable based on the TxBLEND model's daily average salinity results for two example years (2004 and 2005) at the reference node B, a location shown previously on Fig. 1. We employed a spreadsheet approach for calculation of CSD event values with user-selected inputs for the upper and lower range bounds. Although the calculation of CSD could be constrained by seasonal limits, we will use the entire yearly record for this introduction to the concept. The only constraint at the outset is that a string of consecutive days may not continue past the end of a year in order to facilitate the calculations that synthesize multiple year occurrences as presented below. The rationale for, and exceptions to, this end-of-year rule will be explored in the "Methods testing" section.
Salinity duration events, intermediate salinity ranges Figure 6 illustrates the derivation of a set of CSD values but with an intermediate salinity range, specifically a nonzero lower bound, and again without limitation to the time of year, except the same constraint for CSD termination at the end of a year as above. These CSD values are for the salinity range of 2-10 at the same reference node B and for the same time period as used previously. The influence of the lower salinity range bound of 2 is clearly evident with much shorter durations for CSD 2-10 and a greater number of events compared to the lower bound of zero example (CSD 0-10 ) above. This is an expected result since we have simply narrowed the salinity magnitude window from 0-10 to 2-10 and the duration of events in the salinity time series, at the same physical location, that falls with this narrower range would have to be shorter.

Return periods
The variable CSD characterizes salinity by determining the duration of an "event" with a magnitude of salinity in a specific range. Since this variable can be calculated for occurrences within each year, and that set of values further selected for the maximum annual values in each year of the available data, we are also able to develop a measure of the frequency with which certain salinity "events" reoccur through procedures developed in the fields of hydrology and meteorology (Hershfield 1961;Chow et al. 1988).
The variable rCSD is introduced to denote the return period, in years, for specific CSD events as previously defined. This calculated variable is a measure of the frequency with which a certain salinity "event" with salinity in a specific range and of at least a specified benchmark duration, reoccurs over a long period of time. Thus, we will denote this new variable by these two defining characteristics, the salinity range and the event duration: In this step, the set of previously derived maximum annual CSD "events" at each physical location, are further characterized as to those whose durations meet or exceed a qualifying benchmark value, such as 120 d. By determining the number of years between such qualifying events, the calculation method can characterize the frequency of reoccurrence as measured by the period lengths (e.g., once per 5 yr, or once per 10 yr, etc.). As such, it is an analog to the frequency   shown near the top of the figures (e.g., 4 yr between 1987 and 1991 at reference node B). The return period rCSD for CSD 0-10 events of 120 d length is, mathematically, the expected value of the period between the events over a longtime frame, preferably a multidecadal period, but obviously depending upon the underlying salinity data source. The shorter the return period, the more frequently occurring the event is. Herein, we calculate the return period by dividing the total number of years in the record, namely 23, by the number of qualifying events, chosen from the series of maximum events of each year, a technique used in hydrology (see Chow et al. 1988).
Thus, for the example of Fig. 7, for CSD 0-10 events of 120 d length, the return period is 2.56 yr (23/9) or in a more formal mathematical form rCSD 0-10=120days at B52:56 years ð Þ In a more colloquial form, we can say that at reference node B, we would expect salinity to remain in the 0-10 range for 120 consecutive days about one time in each 2.56 yr over a long period of record. For reference node C, the result of this same process yields a return period of 3.29 yr, indicating a slightly less frequent occurrence. This is an expected result for this low salinity range (0-10) at a point further removed from the dominant freshwater source. Table 1 presents a suite of rCSD 0-10 values determined for the seven reference nodes previously shown on Fig. 1 and for several benchmark length values (values shown as infinity indicate that no qualifying events occurred and the return period calculation returns a value of infinity). Through the procedure given above, the return period (rCSD) value at any location, for any particular salinity range, and any benchmark duration can be found.

Contour map development
To portray the new CSD and rCSD variables in a spatial context, the point-by-point values as derived above can be used to generate contoured maps. For the variable CSD, a single year is portrayed by using the maximum annual CSD values of a given year. For the variable rCSD, the derived values are used to produce a map which represent the patterning of salinity over a cumulative period, in this case the entire period of record. In this study, contour maps for the variables CSD and rCSD over the spatial domain of the interlinked Guadalupe and Mission-Aransas Estuary system were derived through spatial interpolation with Geographic Fig. 7. Illustration of the derivation of the return intervals for CSD 0-10 events greater than or equal to 120 d in length at the reference node B (location in Fig. 1). Events that are greater than or equal to 120 d in length are denoted by gray bars. Table 1. Derived return periods for the annual maximum series of CSD 0-10 events at several specific benchmark durations for the reference nodes A though G (locations in Fig. 1 Information System software beginning with the derived values at the 147 "contour" nodes of Fig. 3. The underlying assumption of spatial interpolation is that points in the same vicinity are more similar to one another than they are to those farther away (Tobler 1970). Therefore, any location's value is estimated by weighting the influence of nearer points more than those farther away, a method appropriate for fluid parameters within an estuary system. We used the Geostatistical Analyst extension of ArcGIS 10.2 (ESRI, Redlands, California) to perform interpolation to derive continuous surface contour maps. Since there are a variety of interpolation methods (Lam 1983;Luo et al. 2008), we investigated data patterns and tested various spatial interpolation approaches to identify the most appropriate model prior to beginning analyses. We used a test set of CSD values to assess surface modeling efficacy of the following interpolation methods: inverse distance weighted, local and global polynomials, radial basis functions or spline, and kriging. Those models that produced lower mean absolute error (MAE) and root mean square errors (RMSE), and produced surface models that appeared to reflect environmental patterns reasonably well were deemed as most suitable (Yang et al. 2015). Of these, ordinary kriging using a spherical semivariogram model, with a maximum searching neighborhood of two to seven neighbors (k), in four diagonal sectors produced the best results and was used for the creation of all subsequent contour maps of CSD or rCSD using all 147 contour nodes.
Kriging is analogous to regression analysis in that it produces a statistical model to predict data patterns (a spatial surface in this case) and the differences between model-predicted values and observed values are used to assess error. In comparing different interpolation models using the kriging method for any given CSD or rCSD, those models that produced the lowest RMSE, mean error, and MAE and had mean standardized errors (MSEs) near zero were deemed as superior (Yang et al. 2015). Since RMSE is scale dependent, there is no standardized reference for a "good" fit; it can only be used to compare among interpolation models based on the same data scale and sample size. Each contour map developed with this approach was further assessed for accuracy with the validation nodes as discussed in the "Assessment" section.

Salinity duration event maps
The results of the contouring procedure for one suite of CSD values are shown in Fig. 8, a map of the maximum annual CSD 0-10 for the year 2005, which had total annual inflow near the median (2.91 billion m 3 [2.36 million acrefeet], 12 th rank) in the TxBLEND model's 23 yr record covering 1987-2009. Figure 9, in two panels, shows the same variable of maximum annual CSD 0-10 for both the high-inflow year 2004 (6.78 billion m 3 [5.50 million acre-feet], 3 rd rank of 23) and low-inflow year of 2008 (1.04 billion m 3 [0.84 million acre-feet], 20 th rank of 23).
As expected, there are large differences in the maximum duration values of this low-range salinity event (0-10) among these years of highly differing inflow. For instance, at the location of reference node B, the high inflows of 2004 (Fig. 9, left) were sufficient to maintain salinity in this very low range, below 10, for approximately 150 d, while this drops to about 130 d during the median inflow year (Fig. 8).
In the low inflow year of 2008, salinity was only in this low range for a maximum of approximately 6 d at reference node B (Fig. 9, right). In this portion of the estuary system, inflow variability obviously leads to significant variation in period lengths of continuous low salinity. If one examines reference node E, a location near the seawater exchange with the Gulf of Mexico, a significant contrast is seen to location B. At location E, even in a high inflow year the maximum continuous period duration below 10 was only approximately 5 d, unlikely to represent a biologically or ecologically significant change from the low or median inflow years.
While the yearly results are highly variable, one common trend is evident in the spatial behavior of the CSD 0-10 variable despite the great variability in freshwater inflows: the peak CSD values occur near a freshwater source and decrease steadily as one progresses toward the seawater exchange points with the Gulf of Mexico (northeast of reference node G and south of E). As will be shown in the following subsection, this unidirectional behavior is uniquely linked to the zero lower bound of the salinity range utilized here.

Salinity duration event maps, intermediate salinity ranges
Analogous to the spatial portrayal of CSD 0-10 presented above, Fig. 10, in two panels, illustrates the map view of the maximum annual CSD 2-10 values for the median-inflow (2005) and high-inflow (2004) years, respectively. These intermediate salinity ranges, in this case a nonzero lower bound of 2, and an upper limit well below the maximum salinity values in the system, display particularly relevant features in these spatial portrayals. In contrast to the spatial behavior of the CSD 0-10 variable, where the peak CSD values occur near a freshwater source, there is a distinct "bullseye effect" for CSD 2-10 due to the nonzero lower bound, which shifts the peak values of CSD to mid-estuary regions. For example, in the left panel of Fig. 10, portraying CSD 2-10 for the median year 2005, the values peak in the > 120 d bracket in two regions near reference nodes C and F, before declining again progressing into the core of the Guadalupe Estuary or further up into the Mission-Aransas estuary toward the respective freshwater sources. This effect is the result of these "fresher" areas of the estuary system being more likely to fall below the lower range bound of 2 due to freshwater inflows at the river mouths. A similar effect, although with differing CSD 2-10 lengths, is evident for the high-inflow year of 2004 as shown in right panel of Fig. 10. In fact, for the high-inflow year, the duration of salinity in the 2-10 range is as low in the very upper portion of the Guadalupe Estuary as it is near the outlet to the Gulf of Mexico (near point E) although for opposite reasons. The point E near the Gulf outlet, is rarely in the 2-10 range since salinity nearly always exceeds the upper bound.
The mid-estuary maximum results for CSD with intermediate salinity ranges (nonzero lower range bounds and upper range bounds somewhere below the maximum salinity in the estuary) were always the case for all analyses and mapping in this study of the Guadalupe and Mission-Aransas system. It is also highly likely that this will be the result for any estuary system. In these cases, upper estuary areas closer to freshwater sources will always be more likely to fall below the lower nonzero range bound, while in areas nearer the seawater inlets, there is a higher likelihood of salinity rising above the upper limit. This is in contrast to the unidirectional decrease in CSD under all conditions associated with the zero lower bound case as depicted in Figs. 8,9. It appears that the case of CSD values determined with a zero lower range bound are a singular specialized type of result: the peak CSD values always occur near a freshwater source and decrease steadily as one progresses toward the seawater exchange points. For many estuaries, there may be a "symmetrical" and opposite case that would result from utilizing high range bounds (e.g., 30-40). In many instances, the expectation here would be that the highest CSD durations will be found at the ocean or gulf inlets and these will diminish in a unidirectional fashion toward the freshwater sources. However, for many estuaries, there is a very important exception to this potential symmetry: hyper-salinity (> 35) brought about by high evaporation and inflow deficits during drought. Under these conditions, mid-estuary regions again will display maximum CSD values for even high ranges such as 35-40. So, in summary, it appears that the lower range limit of zero does bring about certain and unique CSD results that are always and unequivocally true: CSD values will be highest near freshwater sources and display a unidirectional reduction toward seawater exchange points. While there may often be an opposite or "symmetrical" case for CSD values with high range limits, we do not believe there is a general and completely predictable outcome at this end of the spectrum due to the complications related to mid-estuary hyper-salinity.
Another notable feature of the maps for CSD 2-10 presented in Fig. 10 is that the maximum values of CSD for this very low salinity range are found in the median-inflow year not the high-inflow year. More specifically, in the median year most of the Guadalupe Estuary exhibited maximum annual CSD 2-10 values in the range of 60-120 d whereas for the high-inflow year this duration was only in the upper portion with much of the lower half of the bay in the 30-60 d range. Again these are the net effects of the intermediate range bounds wherein salinity in many areas during the high inflow year falls below the lower salinity limit for some portion of the year.

Return period maps
While the maps above for CSD are developed for individual years, the maps of rCSD depict return periods, and thus

Johns and Heger
Salinity duration and frequency patterns present long-term salinity behavior over a multiyear period, 23 yr in this case. Figure 11 presents the mapped results for rCSD at the 45 d duration for CSD 0-10 . Figure 12, in two panels, presents the contoured results for 15 d and 120 d duration, respectively. Some interpretation of these return period maps and the information they convey is in order since this type of depiction of salinity events and reoccurring salinity patterns may not be intuitive at first. Figure 11 shows that at reference node A, salinity can be expected to remain in the 0-10 range for at least 45 consecutive days quite frequently: about once per 1.25 yr over the long term. By contrast, at reference node G, salinity can be expected to remain in the 0-10 range for at least 45 consecutive days only once per 15 yr (very infrequently). More broadly, large swaths of the lower Guadalupe Estuary and the transitional area approaching the Mission-Aransas Estuary (reference nodes B and C), would be expected to have salinity in the 0-10 range for at least 45 consecutive days fairly frequently over the long term since this area is covered by the > 1.5-2.5 yr return period zone. By contrast, the lower most portion of Aransas Bay (encompassing reference node E), approaching a Gulf of Mexico connection, would be expected to maintain salinity in this low range for a consecutive 45 d period very infrequently. The greater than 20 yr contour represents areas with either one occurrence or none in the 23 period of record used herein. Collectively, Figs. 11, 12 allow for comparison of the effect of event duration on the frequency of reoccurrence. These three maps all depict salinity events of the same magnitude, namely CSD 0-10 , but at varying durations. At any given point in the estuary, the frequency of reoccurrence falls as the event duration increases. This is represented as higher return periods at any given point as the duration increases. In more colloquial terms, it is "easier" for salinity to be maintained in any given range over shorter durations, hence it will occur more frequently.
Another important result is evident in these figures: rCSD 0-10 values, at any duration, depict a unidirectional decrease in frequency as one moves from the freshwater sources toward the Gulf of Mexico connections. In plain

Johns and Heger
Salinity duration and frequency patterns language, the frequency of maintaining salinity in this particular low range, that includes zero as the minimum limit, decreases, for any duration, as one moves away from the freshwater sources. This is also consistent with the individual year trends of CSD 0-10 as presented previously (see Figs. 8,9). As shown earlier, CSD behavior is much more complex for events with a nonzero lower limit to the salinity range, exhibiting features such as mid-estuary maximums. Similarly, the associated return period variable, rCSD can exhibit complex and potentially insightful spatial patterns for nonzero lower range limits. Figure 13, in two panels, illustrates some of these complex spatial patterns for rCSD 10-20 events of duration 120 d and 45 d, respectively. Fundamentally, it is once again clear that as the duration limit is shortened, the frequency of occurrence for events increases throughout the spatial domain (although not as evident in some locales due to the contour intervals). For example, at reference node B in the main body of the Guadalupe Estuary, a 120-d CSD 10-20 event only occurs about once per 5 yr, while a 45-d event occurs about once per 1.25 yr. Also, the nonzero lower range limit results in a spatial pattern with a distinctive "bullseye effect," analogous to that for the individual year CSD maps with nonzero lower range limit (Fig. 10). For example, the 45-d duration event with salinities in the range of 10-20, had the greatest frequency of reoccurrence, at a return period of 1-1.5 yr, over a broad central spatial portion of the estuary system. From there, the return periods increase (frequency decreases) as one moves toward either the predominant freshwater source, the Guadalupe River, or toward the points of salinity intrusion from the Gulf of Mexico.
Another interesting feature shown in Fig. 13 is that the higher frequency of reoccurrence for the CSD 10-20 events occurs near the weaker of the freshwater sources, the Mission and Aransas Rivers flowing into Copano Bay (refer to Fig. 1 for bay location). This is because the upper portion of the Guadalupe Estuary is much more likely to have salinity fall below the 10 range limit than is Copano Bay due to the much higher inflows typically occurring from the Guadalupe River compared to discharges of the Mission and Aransas Rivers (see Fig. 4). The effect is especially pronounced for the 120-d duration events (left panel) wherein an area with the highest frequency of reoccurrence, with return periods of 1.5-2.5 yr, only occurs in Copano Bay (surrounding reference node F).

Methods testing
The method introduced above for determining the return period tabulated events by using the maximum annual value of CSD in each year of the full period of record, in this case 23 yr . This type of data series is known as an annual maximum series. A frequent criticism of utilizing this type of series is that by choosing the largest, or in this case longest event in each year, the series will typically neglect significant events if multiple such events occur in the same year. For example, the second longest CSD value in some year could exceed the maximum of another year. As a consequence, the annual maximum approach elevates the relative importance of lesser events only due to the fact that they are the maximum of a particular year. In the method utilized above, it also forces the accumulation of consecutive days to restart at the change of the calendar year, as was the case for the example in Fig. 5, potentially missing significant events. Another frequently used option in hydrology is to rely on a partial duration series in which all events that are above a qualifying level or "threshold" are tabulated. The series is also referred to in the literature as a peaks over threshold series (Meylan et al. 2012). For example, if using the partial duration approach for determination of rCSD 0-10 for 45 d, for reference node B in Fig. 5, the second qualifying event (at 46 d) in May-June 2005 would enter into the calculation, whereas it does not in the annual maximum series (AMS) approach. Additionally, with the partial duration series (PDS), the two CSD events at the change of the years 2004-2005 are combined into one event.
A comparison of the computational results for the AMS and PDS approaches is summarized in Table 2 for the CSD 0-10 case at select reference nodes B, D, and G, a subset that exhibited a wide range in return periods previously with the AMS (as in Table 1). As evident in this example, for very frequent to somewhat frequent events, namely, those with AMS return periods in the range of approximately 1-5 yr, the utilization of the PDS may be important, depending on the goal of the analyses. Some return periods of the partial duration approach are approximately one half of those of the annual maximum approach, in this example. For less frequent events, those with AMS return periods exceeding 10 yr, the results tend to converge. This result is consistent with those in the field of hydrology (Chow et al. 1988;Stedinger et al. 1993).
For either of these approaches, it is necessary to have a series of independent events. As shown in Figs. 5, 6, short excursions of salinity only marginally outside the range of interest can break a seemingly single event into several shorter ones. While this is an issue in both approaches as employed herein, recognizing independent events is a frequently cited constraint on utilization of the PDS since all of the smaller events become germane in the calculation. Some potential means of addressing this "noisy signal" issue are presented in the "Comments and recommendations" section.

Spatial interpolation model assessments
All our final models of CSD and rCSD had low MAE (all less than 0.4), MSE near zero, and had relatively low RMSE for their given data scales (Table 3). These diagnostic statistics indicate the interpolation and contouring models and their associated parameters are reasonable and are predicting data patterns well.
We also employed another evaluation method to assess whether the interpolation model used to construct the contour maps adequately predicted salinity patterns throughout the estuary system. Our approach was to use a set of 15 reserved nodes, referred to as validation nodes (Fig. 3) to assess model efficacy. In the course of calculations, CSD and rCSD Table 2. Comparison of return periods for CSD 0-10 events derived with the AMS and the PDS approaches at select reference nodes (previously located on Fig. 1 Table 3. Diagnostic statistics for all CSD and rCSD models indicate good spatial interpolation model fit. values were derived at each of these nodes but they were set aside and not used in the creation of the contour maps. Once contour maps were created, these validation nodes and their derived CSD or rCSD values were overlain to assess the validity of the predicted patterns. Maps were deemed adequate if 12 or more validation nodes correctly fell within the derived contour intervals (since classification accuracy 80% is standard practice) and any exceptions were only marginally outside of those lines. In very limited cases, the kriging parameter k was modified as needed (decreased to 2 or increased to 7) until the resulting contour map predicted at least 12 of 15 validation nodes accurately. An example map depicting validation efforts is shown in Fig. 14.

Discussion
To our knowledge, the calculations and contour map presentations of CSD events within a year and those for rCSD events across longer, potentially decadal, time frames are novel approaches in estuarine science. However, these spatially explicit procedures are analogs of techniques long employed to describe the long-term spatial behavior of meteorological variables such as precipitation extremes. A good example is the compilation of contour maps of Hershfield (1961) which portray maximum precipitation events accumulated over a certain standard duration (from 1 h to several days) that reoccur with specific frequencies (e.g., 1 per year through 1 per 100 yr).
We believe that this integration of spatial and temporal information regarding salinity behavior will enhance our understanding of salinity-modulated processes and habitat relationships within estuaries over and above the currently available techniques. The utility that we anticipate in the application of the CSD and rCSD approaches can be organized into several time scales.
For a single season or year, the spatial portrayal of CSD should enhance our capacity to examine the short-term (e.g., seasonal) effects of salinity on observed biologic and ecologic processes such as nekton or plankton community composition, individual species spawning and recruitment success, disease progression, etc. An example application would include evaluations of the reproduction and recruitment success of Rangia clams (Rangia cuneata), an important mussel in upper areas of Gulf Coast and Atlantic estuaries, in relation to the CSD values of a single year in Fig. 14. An example of the validation exercise performed for each map created. In this case, the two highlighted validation nodes differ slightly from the established contours, but only marginally. Overall 13 of 15 validation nodes agreed with the established contours so this map was deemed to be an accurate representation of the salinity pattern.
Also, calculation and portrayal of CSD events in a spatial context should provide for broad-scale comparisons of salinity behavior among seasons or years with differing freshwater inflow conditions or other physical drivers such as tides or wind (such as we presented in Figs. 8-10). Such inter-year comparisons of CSD could provide new insights into ecological processes such as differing spawning and recruitment outcomes.
A particular subset of short-term examinations of salinity would be those associated with "disturbance" events: those of extreme magnitude and/or duration. Examining CSD of a salinity extreme within a given year or season in conjunction with utilization of the long-term rCSD compilation for context, provides a new vernacular for characterizing these rare events, such as inflow extremes, either high or low, and the related salinity response. Either the reduced salinity associated with a large-inflow event or the elevated values of a drought can be characterized in the return period context (e.g., once in 10 yr occurrence). This should provide a basis for assessing such extreme events for their potential salinityinduced biologic and ecologic ramifications, such as abrupt marsh compositional changes (e.g., Kinney et al. 2014).
Over a longer, multiyear term, examination of reoccurring salinity patterns that may constitute significant physiological stress may offer additional enhancements to the current "Habitat Suitability Paradigm" in widespread use in estuarine species management and restoration. The HSP ranks habitat quality by combining physical attribute scores of estuarine zones with a salinity-based score based on observed patterns, often during abundance peaks over a seasonal time scale (e.g., Clark et al. 2004). While the HSP will find areas of highest suitability, additional insights on long-term suitability could be brought to bear via examination of the likelihood of high-stress or lethal events. The Eastern oyster (Crassostrea virginica), for example, thrives in the intermediate 10-20 salinity, especially during the warmest months. While those salinities are most desirable to support oyster growth, maturation, and reproduction (Cake 1983), it is also well established that low salinity (less than 5), even of short duration ( 20-30 d), can result in over 90% mortality in adults (Loosanoff 1953). On the other end of the spectrum, salinity above 20 for sustained periods in warm weather promotes disease growth and parasite attack, often rampant (Andrews and Ray 1988). In the HSP approach, the procedure typically is focused on determining the geographic area in which 10-20 is the average salinity with a warm-period seasonal constraint (e.g., Cake 1983). This could be complimented by examining the likelihood of undesirable extremes of salinity, both high and low, that may limit oyster survival and/or reproduction in that zone.
For example, examination of reference node B on the long-term average salinity map of Fig. 1 reveals that this location would appear to be ideal for Eastern oysters (C. virginica). That conclusion is reinforced with our results for rCSD 10-20 as shown in Fig. 13 with 45 consecutive days in the 10-20 range occurring almost annually and in the same range for 120 consecutive days about one time in 5 yr. However, if one examines Fig. 11, it is also clear that this location is subject to salinities in the 0-10 range for 45 consecutive days fairly frequently, about once per 2 yr. Further exploration of the duration and frequency in the lethal 0-5 range would be needed to fully explore the suitability of this location. Nonetheless, it is clear that salinity event duration and frequency approaches bring to bear novel information that is not available from the common average salinity portrayal. For species or communities with known salinity tolerances, the salinity event duration and frequency approaches provide a basis for examining extremes that may impinge on long-term habitat suitability.
In a similar fashion, but where salinity tolerances are not well established, the CSD and rCSD approaches provide a basis for determining, to some extent, the salinity characteristics that may limit the spatial range of a particular species or community. As the fields of biogeography and ecology inform, the spatial domain of any species is limited by a combination of physical (extrinsic) factors such as temperature, sunlight, and salinity; biotic (intrinsic) factors such as a species' sensitivity to extremes; and ecological factors such as competition, predation, or the presence or absence of suitable food at the appropriate time (Cox and Moore 2016). With very few exceptions, we do not know these limiting factors with any precision, acting individually or in concert. These knowledge gaps are becoming even more apparent in the emerging field of climate change vulnerability assessments, in which predicting species or community sensitivity and adaptive capacity is hindered by insufficient knowledge of what variable(s) constitute current limiting factor(s) in the existing geographic range limits (e.g., Pacifici et al. 2015). For estuaries, both the CSD and rCSD variables may provide a spatial context for examining the role of salinity-based events in sustaining or limiting species and community ranges.
The rCSD variable, because it is a time-integrator, may prove useful for the development of "correlative" biogeographic models relating a species' range to a multiplicity of abiotic variables (Peterson et al. 2011). Whatever the computational framework, we suggest analyses that balance consideration of: (1) frequently favorable salinity conditions that support growth, food availability (for faunal analyses), reproduction, and the avoidance of disease and parasites, and; (2) the frequency and severity of extreme salinity conditions that would prove intolerable, lethal, or otherwise nonsupportive.
For example, over multiyear to decadal time scales, the rCSD approach could provide insights as to the salinity conditions that may support or limit physical habitat zonation across the spatial domain, such as emergent shoreline or submerged vegetation communities of specific composition. The compositional mosaics of emergent wetlands are generally viewed as products of the complex interactions of elevation, inundation, freshwater availability, tidal flushing, and salinity exposure in the soil pore water (Mitsch and Gosselink 2015). Investigating these mosaics as functions of favorable salinities to support specific marsh plant growth and reproduction vs. the countervailing limits that may be imposed by extremes would be amenable to analyses with our CSD and rCSD approaches.
The consecutive day approach and return period calculations and associated spatial portrayal techniques could potentially be of utility in analyses related to other aquatic stressors. There are analogs of cumulative stress related to periods of low dissolved oxygen (Best et al. 2007;Schmidt et al. 2017) and elevated turbidity (Jalon-Rojas et al. 2015) that could be amendable to the approaches outlined herein depending, of course, upon data availability.

Comments and recommendations
These spatially explicit portrayal techniques for salinity behavior based on Consecutive Salinity Day events and return periods of Consecutive Salinity Day events are applicable to any waterbody with an underlying salinity data set of sufficient temporal length and spatial coverage. Thus, some consideration of these time and spatial constraints is in order.
For limited spatial data, a subset of the presented approach would be to derive just CSD and/or rCSD values at any single point or limited set of points for which a suitable data record exists. The derivation of CSD values for salinity ranges of interest could be based on just a few years of data, for instance to support comparisons of salinity characteristics among years with varying hydro-climatic conditions. However, the derivation of return period information requires a longer and continuous salinity record, probably on the order of two decades at a minimum in order to have any confidence in the patterns of more infrequent events (e.g., 1 per 5 yr or less frequent).
The capacity to fully develop the mapped spatial portrayal of CSD and/or rCSD characteristics as performed herein requires a much more spatially detailed set of data to support the contouring exercise. This step will undoubtedly rely on synthetic data to cover the spatial domain. This could be either simulated data from a hydrodynamic and salinity mathematical model, such as the TxBLEND model used herein , or statistically derived values.
While such synthetic salinity data, whether from a mechanistic or statistical approach, is critical for developing the requisite spatial record, such techniques could also be key to developing a longer temporal record to support return period calculations. Even in the case in which a calibrated and validated model such as TxBLEND is available, the return period calculations herein would benefit greatly if the period of simulation could be extended. Primary drivers of salinity behavior such as freshwater inflows, evaporation, and tidal exchange may have longer records that could be utilized to support statistically based derivations of salinity behavior for extended periods in lieu of complete model extension. These period of record considerations are important in this study and in many potential applications of this technique, in which a short overall period of record for salinity, less than three decades, may be available. In the fields of hydrology (e.g., Stedinger et al. 1993) and meteorology (e.g., Hershfield 1961), return period calculations are often based on periods of record for underlying data generally on the order of 50-100 yr. Record lengths of 2-3 decades may be sufficient for deriving frequent to moderately frequent events (1-5 yr return period) with confidence, but this confidence falls with the rarer events, as was the case in the examples used herein, when only 1-2 events occurred in the analysis over 23 yr.
In addition to the underlying synthetic record development for salinity, there may be some benefits in "pretreating" these data before derivation of the CSD variable. As evident in Figs. 5, 6, the very detailed daily salinity data utilized herein exhibit a somewhat "noisy" signal over the course of a few days to a week period due to the changing and competing influences of primary drivers such as freshwater inflow, tides, and winds. While this detailed picture of salinity behavior is perhaps the most comprehensive, it can lead to interruption of runs of CSD such as that evident in April 2005 in Fig. 5. A very minor excursion above 10 at this location broke a consecutive salinity day period of 124 d that may have otherwise extended to something on the order of 170 d. A couple of options for dealing with this would be easy to implement but were beyond this current effort. First would be an adjustment of the CSD calculation to allow for such minor excursions so long as they were bounded in magnitude and duration. Another approach would be to simply apply a smoothing algorithm to the salinity data, such as the 3-d or 5-d running average. Any such adjustment approach should be rigorously assessed for the sensitivity of the overall CSD and or rCSD results.

Aquatic system considerations in spatial interpolation
Spatial interpolation methods are suitable for continuous surfaces with inherent spatial autocorrelation but applications to aquatic systems have unique characteristics that need to be considered. Because aquatic areas are often interrupted by land features such as islands, peninsulas, inlets, etc., a common course of action is to model the surface and then crop out land areas at the end. Even so, two types of problems can arise, edge effects and disproportionate influence from locations that are, in fact, isolated from the area of interpolation with regard to underlying physical processes. Edge effects can occur because there are fewer inputs from neighbors at edges, and this decreased information results in less accurate predictions at these locations (Conolly and Lake 2006). A rough rule of thumb is to consider the outer 10% of the interpolated surface as suspect (Conolly and Lake 2006).
The disproportionate influence issue arises when points that are distant in an absolute sense exert undue influence on the interpolation results. The normal nearest-neighbor weighting should theoretically minimize influence of points if they are on the other side of a sufficiently sized land barrier. Nonetheless, this effect is seen on the west side of Copano Bay in Fig. 11 where edge effects are compounded by disproportionate influence from locations in Aransas Bay. In this analysis of rCSD 0-10/45days , there were 13 sample nodes to the west of reference node F ( Fig. 3) with actual values ranging from 1.278 to 1.643 and a trend toward decreasing rCSD from east to west. However, the resultant interpolation model produced a contour map with increasing rCSD in the edge areas at the west end of Copano Bay. Apparently, a combination of the lower number of neighbor inputs at this edge along with the influence of high return period values from distant neighbors within Aransas Bay to the south influenced the result, even though those points are across a sizeable land barrier. In retrospect, perhaps additional salinity model points in the west end of Copano Bay may have improved model efficacy. Another option is to use an interpolation method that incorporates barriers a priori. We tested this for the rCSD 0-10/45days analysis (kernel smoothing interpolation for barriers) and it produced a more realistic pattern, but RMSE increased (RMSE 5 8.2703) indicating model fit may not be better. Although this problem did not affect most of this study's analyses, it is an issue that needs further investigation.