I-15_Durham Trend Analysis SOP Final 3-30-15 signed

F0 A
-ki M* rimm.T."
s • �
M
CITY OF MEDICINE
CITY OF DURHAM
WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
City of Durham
Stormwater & GIS Services Division
Water Quality Unit
March 30, 2015
City of Durham Stormwater &GIS Services—[phone]
www.DurhamNC.gov/stormw ter
Design/Plan Review—Drainage/Flooding Concerns—F/oodplain Information
Stormwater Public Education —Surface Water Quality

Approval Sheet
Approved by: Date: � �&
oh Cox, Wa er Quality Manager
�j7
Approved by: A(, w Date: L J
Paul Wiebke, Assistant Director for Stormwater & GIS Services
This standard operating procedure was developed by Tetra Tech (1 Park Drive, Suite 200,
Research Triangle Park, NC 27709, 919-485-8278), March 2014.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE H

TABLE OF CONTENTS
APPROVALSHEET............................................................................................................................................... II
1.0 PURPOSE............................................................................................................................................... 1
1.1 CITY OBJECTIVES...................................................................................................................................................1
1.2 SOP DESCRIPTION.................................................................................................................................................1
2.0 POLICY STATEMENT............................................................................................................................... 3
3.0 BACKGROUND ON THE STATISTICS OF TREND ANALYSIS......................................................................... 5
3.1 TYPES OF TRENDS..................................................................................................................................................5
3.2 KEY STATISTICAL ISSUES FOR TIME SERIES TREND ANALYSIS...........................................................................................8
4.0 STATISTICAL METHODS FOR TREND ANALYSIS.......................................................................................11
4.1
LOWESS..........................................................................................................................................................11
4.2
MANN-KENDALL TEST AND SEASONAL KENDALL TEST .................................................................................................11
4.3
REGIONAL KENDALL TEST ......................................................................................................................................12
4.4
WILCOXONRANK-SUM TEST.................................................................................................................................13
4.5
ARIMA.............................................................................................................................................................13
5.0
SELECTED STATISTICAL PACKAGES AND STATISTICAL METHODS............................................................23
5.1 USGS KENDALL PROGRAM...................................................................................................................................23
5.2 JMPTM..............................................................................................................................................................23
6.0 PROCEDURES........................................................................................................................................25
6.1 PRELIMINARY DATA SCREENING.............................................................................................................................25
6.2 SELECTION OF APPROPRIATE STATISTICAL TESTS........................................................................................................28
6.3 GENERAL INSTRUCTION FOR DOCUMENTING ASSUMPTIONS AND DECISIONS...................................................................29
6.4 MANN-KENDALL AND SEASONAL KENDALL TESTS......................................................................................................29
6.4.1 Preparation of Input Files.......................................................................................................................29
6.4.2 Example: TP at City of Durham Station EL1.9EC....................................................................................31
6.4.3 Example: Water Quality Index at EN8.9ER............................................................................................37
6.4.4 Example: Dissolved Oxygen at City of Durham Station NH3.3SC...........................................................41
6.5 REGIONAL KENDALL TEST ......................................................................................................................................44
6.6 LOWESS..........................................................................................................................................................47
6.7 WILCOXON RANK SUM TEST..................................................................................................................................50
6.8 ARIMAANALYSIS...............................................................................................................................................56
7.0 QUALITY CONTROL...............................................................................................................................63
8.0 REFERENCES.........................................................................................................................................65
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 111

(This page is intentionally left blank)
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE IV

1.0 Purpose
1.1 City Objectives
The overarching objective for this Standard Operating Procedure (SOP) is to guide temporal
trend analysis with regard to the City of Durham's water quality monitoring data to support both
regulatory reporting and public outreach. Trend analysis is defined as the use of numeric and
statistical techniques to determine, within a specified level of statistical confidence, whether a
change has occurred over time and to estimate the rate or magnitude of that change. Trend
analysis can be performed for both concentrations and loads for water quality pollutants of
concern such as nutrients (nitrogen and phosphorus species), toxics, sediment, and pathogens
or physical properties of water such as temperature, dissolved oxygen, pH, chlorophyll a, and
turbidity.
The Division of Stormwater and GIS Services is responsible for reporting to the Durham City
Council whether water quality is improving in response to actions being funded and performed
by the City. This SOP is intended to help City staff conduct a variety of trend analyses,
understand the scientific basis behind each analysis, and to reliably report trends in water
quality and pollutant loading for the City. In addition to communicating results to the Council, it
is anticipated that results will also potentially be used to support public outreach materials, MS4
permit compliance documentation, and other reporting associated with regulatory
requirements or City water quality management strategy implementation.
1.2 SOP Description
This SOP is written at a level for someone with a basic understanding of statistics and water
quality data. It is supported by references from cited literature. The SOP begins with background
on the statistics of trends analysis, including methods for detecting both seasonal and long-term
trends. Statistical concepts are discussed using simple explanations. In-depth discussions are
available in the cited sources. The statistical software and statistical methods covered by the
SOP are identified and briefly described, including an overview of which methods are
appropriate for use by the City to address specific questions.
Prior to applying statistical analysis, it is important to review data for possible limitations. The
SOP discusses considerations for screening monitoring data — what to look for regarding data
quality and quantity, and how to assess what data are appropriate and sufficient for supporting
trend. The SOP then provides detailed examples demonstrating the use of the JMP statistical
software (the commercial statistical software package developed by SAS that the City is licensed
to use) and methods using monitoring data from the City and other sources in Durham County.
The examples lay out the specific steps City staff should follow when performing trend analysis.
When analyses are conducted, all assumptions and decisions regarding data and selected tests
should be documented thoroughly.
Data from several monitoring stations are discussed throughout this SOP. Table 1 provides the
monitoring station site IDs and names for reference. Data were accessed in late December 2013,
and were typically available through October, November, or December of 2013.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 1

Table 1. Monitoring Station IDs and Names
Station ID
Station Name
B3670000
Northeast Creek at O'Kelly Church Rd.
EL1.9EC
Ellerbe Creek at Glenn Rd.
EL5.OEC
Ellerbe Creek at Club Blvd.
EL5.6EC
Ellerbe Creek at Midland Terrace between Club Blvd. & Camden Ave.
EL7.9EC
Ellerbe Creek at W. Murray Ave.
EL8.2EC
Ellerbe Creek at Stadium Dr. west of Duke St.
EL10.7EC
Ellerbe Creek off Sprunt St. upstream of tributary on main stem
EL5.5GC
Goose Creek at Camden Ave.
EL8.1GC
Goose Creek at Holloway St. downstream side of culvert
EL7.ISEC
South Ellerbe Creek at Glendale Ave.
EL8.5SEC
South Ellerbe Creek at Onslow St. and Club Blvd.
EL7.6SECT
South Ellerbe Creek Tributary at Knox St.
EL8.6SECUT
UT South Ellerbe Creek at Corner of Foster and Hunt
EN8.9ER
Eno River at West Point On the Eno at Roxboro Rd.
NE2.2NP
Northeast Creek North Prong at Meridian Parkway upstream
NH3.ONHC
New Hope Creek at Erwin Rd.
TFO.OTC
Third Fork Creek at Hwy 751.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 2

2.0 Policy Statement
The technical procedures of this document are used by staff that have training in statistics and
that have experience in advanced statistical analysis.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 3

(This page is intentionally left blank)
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 4

3.0 Background on the Statistics of Trend Analysis
3.1 Types of Trends
Time series data may exhibit a variety of patterns, some of which contain systematic trends.
Data may exhibit random variability or seasonal patterns without any systematic trend, or there
may be consistent trends with time imposed on the underlying variability and/or seasonal
cycles. Time series behavior can also change suddenly, signifying a change in the environment
affecting the parameter. Several types of time series behavior and the circumstances causing
the behavior are discussed below. Graphs depicting the types of trends and other time series
patterns are shown in Figure 1.
Nn Trend
Time series data may exhibit no systematic trend and represent a series of measurements that
fluctuate randomly around a constant mean. A key challenge for time series analysis is
distinguishing the absence of trend from a trend that is obscured by random noise. This is more
likely to be an issue when data from a relatively short time period are evaluated; trends may be
gradual and cannot be discerned from random stochastic processes without more observations.
For water quality, the absence of a trend is most likely when the conditions causing pollutant
loading do not change over time except for the variability associated with weather patterns
(assuming that weather itself is not subject to long-term trends).
Monotonic Trend
The most obvious time series trend is a monotonic trend. A monotonic trend is one that
represents a change in a single direction over time, often superimposed on other types of
variability. Monotonic trends may be increasing or decreasing, and may be linear or nonlinear.
For example, land use change over time often influences the load and concentrations of
pollutants. Management measures implemented gradually over a period of time can result in
monotonic trends, such as stormwater infrastructure upgrades and stringent stormwater
management requirements applied to re -development in an older urban watershed. Even
remote actions can influence water quality; for instance, nitrate wet deposition in the Durham
area appears to be decreasing since the early 2000's (Tetra Tech, 2013), likely reflecting national
efforts to control oxidized nitrogen emissions from coal-fired power plants.
Seasonal Pattern
The concentrations of many pollutants vary systematically throughout the year and it is
important to distinguish these seasonal patterns from systematic multi -year trends. A seasonal
pattern varies in a cyclical fashion over the course of a year. Seasonal patterns may be quite
obvious and expected — such as the one shown in Figure 2, which represents dissolved oxygen
(DO). DO is highly influenced by temperature and thus has strong seasonality, with the highest
values occurring in the winter and lowest values in the summer. Sometimes seasonal patterns
are present but less obvious; for instance, nitrate concentration often increases during the
springtime in locations where fertilizer application is prevalent.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 5

No Trend
10
9
8
7
6
5
4
3
2
1
0
Jan Apr Jul Oct Jan Apr Jul Oct Jan Apr Jul Oct
10
9
8
7
6
5
4
3
2
1
0
Jan Apr Jul Oct Jan Apr Jul Oct Jan Apr Jul Oct
Monotonic Trend
10
9
8
7
6
5
4
3
2
1
0
Jan Apr Jul Oct Jan Apr Jul Oct Jan Apr Jul Oct
Seasonal Pattern I I Seasonal with Trend
Step Change
10
9
8
7
6
5
4
3
2
1
0
Jan Apr Jul Oct Jan Apr Jul Oct Jan Apr Jul Oct
Trend Start
10
9
8
7
6
5
4
3
2
1
0
Jan Apr Jul Oct Jan Apr Jul Oct Jan Apr Jul Oct
10
9
8
7
6
5
4
3
2
1
0
Jan Apr Jul Oct Jan Apr Jul Oct Jan Apr Jul Oct
Impulse
10
9
8
7
6
5
4
3
2
1
0
Jan Apr Jul Oct Jan Apr Jul Oct Jan Apr Jul Oct
Serial Correlation
10
9
8
7
6
5
4
3
2
1
0
Jan Apr Jul Oct Jan Apr Jul Oct Jan Apr Jul Oct
Figure 1. Examples of Time Series Patterns (adapted from Gilbert, 1987)
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 6

Figure 2. Time Series of Dissolved Oxygen at City of Durham Station EL7.9EC
Seasonal Pattern with Monotonic Trend
Seasonal patterns and long-term trends are often superimposed, which complicates trend
detection. DO could have an underlying monotonic trend if there were a gradual change in
upstream influences. For instance, biochemical oxygen demand (BOD) loads might increase
downstream of a WWTP over time, leading to changes in DO in addition to the expected
seasonal pattern. Watershed development could also influence DO, stemming from changes in
baseflow and a flashier flow regime.
Step Change
A step change occurs when there is a sudden shift in the magnitude of a pollutant. This is often
seen downstream of a WWTP following an upgrade in treatment technology. Phosphate
detergent bans have resulted in rapid declines in observed phosphate and TP. A change in
laboratory methods (or a change in the laboratory used to perform the analysis) could also lead
to a step change in the reported data, even if the true values of the parameter have not shifted.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 7

Impulse
An impulse represents a temporary shift in the magnitude of a pollutant. This can occur when a
new, but temporary pollutant source with a strong signal occurs upstream. For example,
construction activity may elevate TSS and turbidity readings until the project is completed and
the site stabilized. A leaking sewer line could increase downstream bacteria levels until the leak
is detected and the repair made. Distinguishing a temporary impulse from a monotonic trend or
a step change can be challenging, especially in the period prior to the cessation of the impulse.
Beginning of a New Trend (Trend Start)
As noted previously, a monotonic trend represents a change in a single direction over time. If
there is a new trend (or a change in trend), identifying the point of change becomes the
question of interest. This is important to point out because different statistical tests are needed
to detect changes or shifts in a trend (i.e., some form of inflection point). In the case of the
example graph, the inflection point is clear and one might simply select an arbitrary point to
separate the data into two series for analysis. However, in other cases the inflection point may
not be clear — or one may need to test that there is a change to begin with. Graphical screening
and use of outside knowledge are a critical component of the selection of the proper statistical
test, and are discussed further in Section 6.1.
Serial Correlation/Autocorrelation
Classic regression models are founded on several underlying assumptions about the data: that
the data are normally distributed, that the errors have a constant variance, and that deviations
from the mean or from the underlying trend in consecutive data points are independent of each
other. Time series monitoring data often break the last assumption — a particular measurement
may have a similar magnitude to adjacent observations in time. This is referred to as serial
correlation or autocorrelation. Monitoring data often show positive serial correlation, meaning
that positive deviations (about the mean or trend) at one point in time tend to be associated
with positive deviations in adjacent time periods (and negative errors are associated with
adjacent negative errors). The example graph (Figure 1) shows an example of positive serial
correlation. An important challenge for time series analysis is distinguishing serial correlation
from true underlying trends in concentration. The problem of serial correlation is discussed in
the following section.
3.2 Key Statistical Issues for Time Series Trend Analysis
Analysis for trends may use either "parametric" tests (tests that require assumption of a specific
form of the statistical distribution) or "non -parametric" tests (which do not require the
assumption of the form of the distribution). In discussions of trend analysis, the first method
commonly covered by statistical texts is linear regression, which is a parametric method based
on the assumption that the data follow a normal distribution. On the surface, this would seem
an ideal way to test for trend, with time as the independent variable and concentration or load
as the dependent variable. However, pollutant monitoring data typically exhibit properties that
confound the use of linear regression. More often than not, data violate two of the core
assumption of linear regression: the data are not normally distributed, and the data are serially
correlated. As a result, at a given confidence level, it is possible to have a finding of a trend (i.e.,
reject the null hypothesis of no trend) when in fact no trend is present. Linear regression could
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 8

be used if the data have a normal distribution and do not have serial correlation. Generally
speaking, linear regression is more efficient (i.e., more likely to detect a measure of statistical
significance) than non -parametric procedures if the proper assumptions are met. However, non -
parametric procedures have much higher efficiency when data have even small departures from
normality (Hirsch et al., 1991), which is the case more often than not for water quality
monitoring data. In addition, a detection limit that is close to (or higher than) the ambient
magnitude of the pollutant presents a further challenge for using linear regression. For these
reasons, linear regression is not considered a useful method for time series trend analysis in this
SOP.
Seasonal variation in data is also commonly encountered, and may be difficult to discern in
graphical plots of data. The pattern is often unmistakable in DO data, but other parameters may
exhibit variation with less of a magnitude shift. The shift may also occur during only a few
months of the year. Seasonality complicates trend analysis by making it more difficult to detect
underlying trends. This problem can be addressed by using statistical methods that remove the
seasonal cycle from the data prior to trend analysis, or methods that simply account for
seasonality as part of the overall test.
Pollutant concentrations are often correlated to the rate of streamflow. This occurs primarily for
two reasons: dilution, and storm event washoff processes. If a pollutant is delivered to a stream
at a relatively constant rate, such as with a point source or as part of groundwater, then a higher
volume of flow will dilute the pollutant. On the other hand, during storm events a pollutant may
be washed off land surfaces in higher concentrations than the background concentration in the
receiving stream. Wet weather collection system discharges can also increase concentration
during high intensity storm events. The result after mixing may be an increase in observed
concentration with high flows. A statistical method for addressing variance in pollutant data due
to differences in flow is discussed in the following section.
Another key issue in trend analysis is distinguishing random variability from true trends. It is
desirable to avoid uncertain claims of trend that have a risk of being proven untrue as further
data are collected. The analyst should therefore conclude that a trend is statistically significant
only when the estimated probability or confidence level that the trend coefficient is different
from zero is 95 percent or greater (i.e., there is less than a 5 percent risk that a trend is
incorrectly identified when none actually exists). The choice of a 95 percent confidence level is a
traditional but essentially arbitrary choice regarding an acceptably low risk of identifying a trend
where none exists. If another confidence level is used, it should be agreed upon and
documented by staff members prior to conducting tests.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 9

(This page is intentionally left blank)
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 10

4.0 Statistical Methods for Trend Analysis
4.1 LOWESS
The efficiency of procedures for detecting trends can be improved whenever extraneous
variance in data is reduced. As noted in the previous section, pollutant concentrations are
frequently correlated with flow. The presence of flow effects on concentration may mask a
trend, or exaggerate it; one cannot tell from a concentration time series whether an apparent
trend is due solely to systematic differences in discharge over time (for instance, a series of wet
years), or truly due to a real trend in the concentration. One method for reducing extraneous
variance is called LOWESS, which stands for locally weighted scatterplot smoothing (Cleveland,
1979). A smoothing window with a user -defined width in terms of the fraction of total
observations is passed over the data, one observation at a time. The smoothing function gives
stronger weight to observations with less variance compared to the center point, while
observations with greater variation have less weight. As a result, LOWESS is less sensitive to
outliers and is preferable over parametric regression procedures for time series analysis. In the
case of concentration and flow, LOWESS is performed on the relationship between the two
(Hirsch et al., 1991). Once the variation due to flow is removed, trend analysis is conducted on
the residuals. The choice of the width of the smoothing window (denoted as fl may be
important; the goal is to select f to minimize variability of concentration related to flow, while
still allowing any underlying trend to remain undistorted. Cleveland (1979) recommends using a
value of 0.5 for f as a starting point, with variations in the range of 0.2 to 0.8 being reasonable.
4.2 Mann -Kendall Test and Seasonal Kendall Test
The non -parametric Mann -Kendall test for trend (Mann, 1945; Kendall, 1975) forms the basis of
the method that is likely the most frequently used for trend analyses performed on water
quality monitoring data —the Seasonal Kendall Test. The method was developed and
popularized by USGS researchers throughout the 1980s (Hirsch et al., 1991), and USGS published
computer code supporting its use.
Mann -Kendall is especially useful for detecting trends in environmental variables for several
reasons:
• The test is non -parametric, and the data do not need to be normally distributed.
• Missing values are allowed; gaps are simply ignored.
• Data reported at the detection limit can be used without censoring, so long as the
values are set lower than the smallest observation.
• Short periods of record are allowed, though small sample sizes will affect the
significance of the result.
This is all possible because Mann -Kendall looks only at the relative magnitudes of sequential
data, so the type of distribution, gaps, and the assumptions used for non -detects become
irrelevant. The test does, however, assume that the deviations about the trend for the observed
data are not strongly serially correlated. This is likely to be the case for infrequent (e.g.,
biweekly) observations even if daily observations do show serial correlation.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 11

The Mann -Kendall test is performed by first creating an indicator function that takes values of
-1, 0, and 1. The data are ordered in time, and -1 is assigned when the value at time t+1 is less
than the value at t; 0 is assigned when the values are equal, and 1 is assigned when the value at
time t+1 is greater than the value at t. The procedure is repeated at time t+2, then t+3, and so
on until all possible comparisons are made. The Mann statistic S is then calculated as the sum of
the values from the indicator function. If S is a large positive value, then measurements taken
later in time tend to be larger than values taken earlier in time (and the converse is true for a
large negative S). The statistic is compared to tabulated values for smaller datasets; when n
exceeds 40, the variance of S is calculated using a complex formula, a Z statistic is calculated
from S and its variance, then Z is compared to the standard normal distribution. Some statistical
texts discuss (and programs report) Kendall's tau rank correlation coefficient, which is an
alternate measure equivalent to Mann's S statistic, and is used to arrive at the same Z statistic.
The Seasonal Kendall test is a generalization of the Mann -Kendall test that tests for trends
underlying seasonal patterns, and was developed by Hirsch et al. (1982). In its original
application, data were divided into 12 "seasons," with each month representing a season.
Multiple observations within one season are handled by taking the median of the values. The
Mann -Kendall test statistic and its variance are calculated separately on each season. The
statistics are summed and a Z statistic computed. The null hypothesis Ho is there is no trend,
while the alternative hypothesis HA is either an upward or downward trend (a two -tailed test). Z
is compared to tabular values for smaller sample sizes, or the standard normal tables are used
for larger sample sizes. Hirsch and Slack (1984) note the test can be performed on as little as
two years of monthly data, though fewer data will of course affect the significance of the result.
Serial correlation among values within a season can be addressed by a modification of the test
statistic (Hirsch and Slack, 1984). The modification is recommended in cases where there are 10
or more observations per season (i.e., 10 years of data if seasons are defined monthly) due to
difficulties in accurately determining covariance for fewer data.
A rate of change or trend slope can be calculated as well for the Seasonal Kendall test. The slope
is based on Sen's non -parametric slope estimator (Sen, 1968). The method estimates a series of
slopes between values from the same season. The seasonal Kendall slope is the median of this
series of slopes. A confidence interval can also be estimated.
4.3 Regional Kendall Test
The Regional Kendall test (Helsel and Frans, 2006) is a recent development in USGS's continuing
work in trend analysis. It is an extension of the Seasonal Kendall test; instead of data being
divided into different seasons, unique locations are tested in conjunction with each other in
place of seasons. To be meaningful, the stations should be from the same "region", i.e., have a
spatial area in common, but the test does not account for specific location (in other words, it is
not a spatial test that uses coordinates). The null hypothesis Ho is that there is no overall trend
among the stations, while the alternative hypothesis HA is the presence of an overall upwards or
downwards trend among the stations. As is the case with Seasonal Kendall, if some locations
show an upward trend while other show a downward trend, the trends cancel each other out
and the null hypothesis cannot be rejected. The test looks at the overall pattern across
locations; a particular location may have a weak trend, but contribute to a finding of overall
trend across multiple stations. It is important to note that failing to reject the null hypothesis
does not signify that no trend is present at any station. It is possible for no trend to be present
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 12

at any station, but it is also possible that upward and downward trends at different stations
have canceled each other out.
Helsel and Frans recommend using a minimum of five years of monthly observations at each
location to have sufficient data to conduct the test. In addition, very closely spaced locations
should be avoided, because measurements are more likely to be replicates of each other than
independent measures of the process. Multiple values for a given year and location are allowed;
the median is used as is the case for Seasonal Kendall when multiple values are present for a
given year and season.
4.4 Wilcoxon Rank -Sum Test
The Wilcoxon Rank -Sum test, also known as the Mann -Whitney U test and the Wilcoxon-Mann-
Whitney test (Mann and Whitney, 1947), is a non -parametric test of the relationship of central
tendency between two independent (non -paired) data sets of sizes nl and n2. It is similar to the
independent sample t-test, but can be used with non -normal distributions and data sets with a
moderate number of non -detects. The test is useful for detecting a step change when the point
in time the change occurred can be identified (such as a WWTP upgrade or the implementation
of a management measure). The data are simply divided at the change point into two sets, and
the test is performed. The test does not determine the magnitude of the change, but is efficient
for determining whether means differ between data sets with non -normal distributions. Gilbert
(1987) provides a good outline for performing the test.
The null hypothesis (Ho) is that the populations from which the two data sets derive are the
same, whereas the alternative hypothesis (HA) is that the populations have different means. The
procedure is complicated (and the formulas differ depending on the presence of ties), but
entails pooling the data as m = nl + n2, ranking them, and assigning them values of 1 through m.
Next, the ranks of the nl data only are summed. This sum is used in conjunction with nl, n2, and
m to calculate the test statistic. The statistic is compared to tabled values of nl, n2, and a for
small sample sizes, or to the normal distribution for larger samples sizes. There is no lower limit
on the number of data needed to conduct the test, but of course small sample sizes will affect
the significance of the result.
The Hodges -Lehman estimator of step -trend magnitude (Hodges and Lehman, 1963) is often
associated with the Wilcoxon Rank -Sum test. The estimator is calculated as the median of all
possible differences in the data between the two time periods.
4.5 ARIMA
Autoregressive integrated moving average (ARIMA) models provide a powerful toolset for
interpreting and forecasting time series based on a history of observations. ARIMA is built on
the use of stochastic process time series models. The approach was originally published by Box
and Tiao (1975) and later described by Box and Jenkins (1976); sometimes the models are
referred to as Box -Jenkins models. There are (up to) three parts of an ARIMA model — the
autoregressive portion (AR), the moving average portion (MA) and the time series differencing
component (1). ARIMA models may have one, two, or all three elements in various
combinations, and are generally written in the following form:
ARIMA(p,d,q)
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 13

where p is the number of the autoregressive term, d is the number of differences, and q is the
number of the moving average term as discussed below in this section. Models that omit one or
more of the terms are often shown as abbreviations of ARIMA with p, d, or q omitted as
appropriate. For instance,
• AR(1) — Autoregressive of order 1
• MA(2) — Moving average of order 2
• ARMA(1,1) — Autoregressive of order 1 and moving average of order 1 (no time series
differencing)
• ARI(1,2) — Autoregressive of order 1 and time series differencing of order 2
Pindyck and Rubinfeld (1997) provide a solid foundation for understanding ARIMA models and
time series analysis in general, although the focus is from an econometric perspective. Chandler
and Scott (2011) include a chapter on stochastic process models including ARIMA, with a focus
on applications to the environmental sciences. An excellent reference for ARIMA modeling is the
online class notes for an economic forecasting class from Bob Nau of Duke University (Nau,
2005).
In the following sections, aspects of ARIMA modeling are briefly discussed. First, the critical
issue of stationarity of data series is introduced. Next, each component of an ARIMA model is
discussed to provide a general understanding. All the components of ARIMA can be combined in
various forms, and formulae become progressively complicated as components are added and
order increased. The section concludes with the steps that should be undertaken to develop an
ARIMA model.
Stationarity and ARIMA
The autoregressive and moving average portions of ARIMA are calculated on stationary time
series. A time series that is stationary has constant statistical properties over time; mean,
variance, and autocorrelation do not vary over the course of time, even though individual
observations may vary. This property is termed stationarity. If the statistical properties do vary
overtime, the series is nonstationary. Astationarized time series has been transformed in some
manner to create a stationary time series from a nonstationary time series. The "integrated"
part of ARIMA can provide the transformation needed to stationarize a time series with respect
to trends in the mean and variance. Note that transformation of data with a log -normal
distribution (e.g., pollutant loads) may be necessary prior to running the ARIMA analysis in order
for the time series to be stationary in terms of its variance. In addition, ARIMA models are built
on the assumption that data are spaced equally in time (e.g., monthly), and that the time series
is complete. A few missing values can be tolerated, but there should be no major gaps in the
dataset. A minimum of fifty observations should be available to perform ARIMA, though more
data are preferable (Gilbert, 1987).
Autocorrelation in time series is addressed through the autoregressive and/or moving average
parts of the ARIMA method. The autocorrelation function (ACF) and partial autocorrelation
function (PACF) can be used to diagnose presence of autocorrelation in a time series. The
autocorrelation function provides the degree of autocorrelation between data points at t + k
lags, where t is each observation and k is the series of lags 1...n. By default JMP displays the ACF
up to 25 lags, blue lines represents two standard errors (Figure 3). If the autocorrelations for
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 14

k > 0 are all close to zero, the series is likely stationary, which is the case in the example shown
in Figure 3. The partial autocorrelation function (PACF) provides the correlation between a
variable and lags of itself, and signifies the amount of correlation not explained by correlations
at lower lags (Figure 3). Both the ACF and PACF are useful tools for specifying ARIMA models,
and are discussed in more detail at the end of this section. The Ljung-Box Q statistic
communicates accumulated sample autocorrelation up to the reported lag; the p-value is an
indication of whether the time series shows significant autocorrelation across the preceding
lags. The ACF/PACF plots shown in Figure 3 are for example only to show the analysis of
autocorrelation in a stationary time series. Numerous values below the detection limit are
present in the dataset, limiting its utility for further ARIMA analysis.
Time Series NOx
1
Mean
0.2421857
0.8
Std
0.1627216
N
70
O 0.6
Zero Mean ADF
-3.562136
z 04
Single Mean ADF
-7.661944
Trend
ADF
-7.613366
0.2
0
01/01/2007
01/01/2009 01/01/2011
01/01/2013 01/01/2015
Date
Time Series Basic Diagnostics
Lag
AutoCorr-.8-.6-.4-.2 0 .2 .4 .6 .8
Ljung-Box Q
p-Value
Lag
Partial-.8-.6-.4-.2 0 .2 .4
.6 .8
0
1.0000
0
1.0000 ;
1
0.0556
0.2258
0.6347
1
0.0556
2
0.0707
0.5960
0.7423
2
-0.0740
3
-0.0671
0.9348
0.8170
3
-0.0593
4
-0.0798
1.4213
0.8405
4
0.0787
5
0.0202
1.4529
0.9184
5
0.0200
6
-0.0280
1.5147
0.9585
6
-0.0462
7
-0.0429
1.6617
0.9762
7
-0.0466
8
0.0409
1.7974
0.9866
8
-0.0468
9
-0.1669
4.0979
0.9049
9
-0.1749
10
0.0733
4.5491
0.9192
10
0.0737
11
0.0288
4.6199
0.9482
11
-0.0185
12
0.0164
4.6432
0.9688
12
-0.0016
13
0.0720
5.1018
0.9729
13
-0.0993
14
0.0647
5.4789
0.9780
14
0.0923
15
-0.0063
5.4825
0.9872
15
-0.0505
16
-0.0632
5.8555
0.9896
16
-0.0733
17
-0.0648
6.2551
0.9913
17
-0.0785
18
0.1036
7.2947
0.9874
18
0.0991
19
0.0050
7.2972
0.9925
19
-0.0275
20
-0.0130
7.3142
0.9955
20
-0.0241
21
-0.0765
7.9170
0.9955
21
-0.0839
22
-0.0574
8.2635
0.9964
22
-0.0696
23
-0.0465
8.4951
0.9974
23
0.0397
24
0.1350
10.4926
0.9922
24
0.1061
25
0.2431
17.1107
0.8777
25
0.2086
Figure 3. ACF and PACF for NOx at City of Durham Station NH3.3SC
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 15

Autoregressive Models
In an autoregressive model, the value of the variable of interest at any given time is dependent
on the value of the variable at one or more previous time steps. An autoregressive model (of
order p) takes the following form:
Yt = OlYt-1 + 02Yt-2 + ... + Opyt-p + 8 + Et
where yt is the current observation, p is the order, 01 to yip are a series of weighting coefficients,
6 is an optional constant term, and Et is a random disturbance process. yt is generated from the
weighted average of past observations going back p periods, with unique weights 01 to Op for
each period. Stated differently, p is the number of terms in the model that describe dependency
on successive observations. An AR(2) model is predicted by the two previous observations.
Moving Average Models
While an autoregressive model is dependent on a combination of previous observations and a
noise term, moving average models are based on a weighted average of the previous noise
terms. It takes the following form:
yt = µ + Et - O1 Et-1 - 02 Et-2 - ... - OQ Et—Q
where yt is the current observation, q is the order, 01 to Op are a series of weighting coefficients,
µ is an optional constant term, and Et is a random disturbance process. As is the case for
autoregressive models, the random disturbances are assumed to be independently generated —
in other words, normal random variables with a mean of zero and no covariance. Another way
of thinking of a moving average model is that it is dependent on the persistence of the
disturbance, or random shock from one observation to the next. In an MA(2) model, the current
observation depends on the variation of the two previous observations.
Integrated Models (Time Series Differencing)
In time series differencing, a new time series is created as the difference between successive
observations. Time series differencing represents the "integrated" component of ARIMA, and is
used to take a nonstationary ARMA model to make it stationary. An ARIMA model with one
order of differencing (without an AR or MA component) would be written as ARIMA(0,1,0), or
also as I(1). More than one order of differencing can be specified, where a second order
differencing is given by the difference between first order differences. A first -difference
integrated model estimates the stationary series, y, from the untransformed series, x, as:
yt = xt - xt-1 + it + Et
where is a constant term, and Et is a random disturbance process with expected value of zero.
Time series differencing removes trend to make a series stationary; thus the need for an "I"
component in the ARIMA model is indicative of the presence of a deterministic trend that
cannot be explained by the random components introduced by the autoregressive and/or
moving average processes.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 16

Seasonal ARIMA
Each component of an ARIMA model can have a separate seasonal order— i.e., for monthly data
at lag 12. The nomenclature for seasonal ARIMA is as follows:
ARIMA(p,d,q)x(P,D,Q)12
where P, D, and Q are the seasonal orders respectively for the autoregressive, differencing, and
moving average components. Both regular and seasonal orders can be specified in any
combination. The "12" connotes a lag of length twelve, appropriate for monthly data. A
different seasonal lag could be used, but is not addressed (nor needed) in this SOP.
Guidelines for Soecifvina an ARIMA model
There are three steps in the development and application of an ARIMA model:
1. Identification, where various models are tested for stationarity and various orders of
stationary ARIMA models are identified.
2. Estimation, where the candidate models are compared to judge their adequacy.
Forecasting, where the selected models are used to predict future values and/or
estimate trend.
The identification stage is by far the most complex, and requires a great deal of judgment. Most
statistical texts and guidance recommend having a solid understanding of the underlying ARIMA
models and the influence of various types of ordering. The form of the model
prediction/estimation of trend depends heavily on how the model is ordered. The following
steps/rules are paraphrased from Nau (2005). The instructions assume that successive models
will be examined for stationarity. In addition, no more than two orders of any term should ever
be considered. The minimum number of years of data needed to fit an ARIMA model is four to
five.
Differencing (i.e, the d term in ARIMA). Note that the total order of differencing is the
sum of regular plus seasonal differencing (d+ D).
a. The correct amount of differencing is the lowest order that results in a
stationary time series — i.e., one that fluctuates around a well-defined mean,
and whose ACF decays rapidly to zero. If the series has positive autocorrelation
out to a high number of lags, then an order of differencing needs to be
introduced or added.
b. Conversely, if the first lag (k = 1) autocorrelation is zero or negative, or if all the
autocorrelations are small and have no pattern, then a higher order of
differencing is not needed. If the first lag autocorrelation is -0.5 or more, the
series is likely over-differenced, and the order of differencing should be
reduced.
c. The best differencing order is likely the one where the standard deviation is the
lowest.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 17

d. A model with no differencing assumes the original series is stationary (and thus
has no trend other than what is attributable to random processes). One
difference indicates the series has a constant average trend. Two differences
means the series has a time -variable trend.
Include a constant term (i.e, intercept) when specifying models with no
differencing or one order of total differencing. Omit the constant in models with
two orders of total differencing.. For a one -difference model the constant
represents a deterministic linear trend.
2. Identifying the order of the AR (p) and MA (q) terms.
a. If the PACF shows a sharp cutoff while the ACF shows a first lag positive (and
significant, i.e., above two standard errors) autocorrelation, then an AR order
should be considered. The lag where the PACF shows the sharp cutoff is the
number of the AR order.
b. If the ACF shows a sharp cutoff and/or the first lag autocorrelation is
significantly negative, then an MA order should be considered. The lag where
the ACF cuts off is the number of the MA order.
c. If a model has both an AR and an MA order, try a model with one fewer of each
of the terms. It is possible for AR and MA orders to counteract each other.
d. If the sum of the AR coefficients (including their standard errors) reported by
the model is nearly equal to one, then the model is said to have a "unit root." If
this occurs, reduce the AR order by one and increase the differencing by one.
e. If the sum of the MA coefficients (including their standard errors) reported by
the model is nearly equal to one, this also indicates presence of a unit root. If
this occurs, reduce the AR order by one and reduce the differencing by one.
f. If long term forecasts reported by the model are erratic or unstable, there may
be a unit root in the AR or MA orders.
3. Seasonal terms.
If the series and ACF show a strong and consistent seasonal pattern, then use an
order of seasonal differencing (D). Never use more than one order of seasonal
differencing or more than two total orders of differencing (d + D).
b. If the ACF shows an overall tendency towards positive values, adding a seasonal
AR order (P) may help. Conversely, if the ACF shows an overall tendency
towards negative values, adding a seasonal MA order (Q) may help. Do not mix
seasonal AR and MA terms in the same model, and avoid two orders of seasonal
AR or MA within a model.
During the estimation phase, models should be compared to each other to determine the best
candidate for forecasting or trend estimation. JMP provides the Model Comparison window,
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 18

which provides several metrics for each model considered and a graph of the model with
confidence limits. An example of a model comparison window is shown in Figure 4 for NOx at
the City of Durham station EL1.9EC. Generally speaking, models with lower variance are
preferable. The Akaike information criterion (AIC) is a good measure of overall performance; the
lower the AIC, the better the model is likely to be. When measures of model performance are
nearly equal, the simpler model is always preferable to the more complex one.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 19

Report
Graph
Model
DF
Variance
AIC
SBC RSquare
-2LogLN
Weights .2.4.6.8
MAPE
MAE
2
2
- IMA(1, 1)
113
2.0532372
411.76578
417.25565
0.205
407.76578
0.383289
91.504529
0.987610
2
0
- IMA(1, 1) No Constrai
113
2.0532372
411.76578
417.25565
0.205
407.76578
0.383289
91.504529
0.987610
2
❑
- ARMA(1, 1)
114
1.9490964
413.48797
421.77449
0.244
407.48797
0.162017
95.757386
0.9687
8
2
11
- AR(3)
113
1.9820658
416.40398
427.45267
0.238
408.40398
0.037701
95.136516
0.983011
2
0
- AR(2)
114
2.0056766
416.75742
425.04394
0.222
410.75742
0.031594
99.075009
1.0072
0 o
2
0
- AR(1)
115
2.1209753
422.17028
427.69463
0.171
418.17028
0.002110 :
106.72285
1.0632
5 Z
2
0
- I(1)
114
2.9870101
453.19288
455.93781
-0.16
451.19288
0.000000
102.39818
1.1916
3
Figure 4. Example Model Comparison Window in JMP
g
7
6
5
4 .j
3
2 ,�•�• �'�...
0
01/01/2003 01/01/2007 01/01/2011 01/01/2015
1.0 Date 1.0
0.5 0.5
0.0 - - - 0.0
-0.5 -0.5
-1.0 1.0
0 5 10 15 20 0 5 10 15 20
Residual ACF Residual PACF
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 20

In the forecasting stage, the model is evaluated as an instrument for forecasting future
conditions, and/or eliciting trend. The apparent trend predicted by the model is equivalent to
the direction and degree of change in the model forecast —which may include both a
deterministic trend in the mean and a random component due to the persistence of shocks. The
net change from the present period to a forecast future period can be found by exporting the
results of an ARIMA model as follows:
• Select the down —pointing red triangle next to the Model, and choose Save Prediction
Formula.
• A new window appears (Figure 5), which includes predicted values.
• Use the forecasting period (past the observed data time period) to estimate the future
trend by differencing values.
• If the forecast has a seasonal component, use months spread a year apart to find the
difference.
The non-random or deterministic trend in an ARIMA model is given by reversing the time series
differencing (d) portion of the ARIMA model, which is equivalent to restoring the non-
stationarity in the mean. This can be obtained by summing the stationarized time series a total
of d times. Over the entire model period, the deterministic trend is given by the change
between the last value in the undifferenced series (at time t) and the initial undifferenced
observation, where the last value in the undifferenced series is given by:
xt = xp + Zit=i Yi,
where x represents the undifferenced series and y represents the differenced series. Because
each realization of yr is formed from the previous value of x, a random disturbance with
expectation zero, and a constant intercept term, the intercept fit in the model is an estimate of
a linear deterministic trend in an ARIMA(p,1,0) model. (For the more complex cases of d > 1 and
mixed models with both d and q > 1 consult Pindyck and Rubinfeld (1997) or other advanced
time series modeling texts.) If the ARIMA model does not require differencing (i.e., it is an ARMA
model with d = 0) then the model does not contain a deterministic trend.
NOTE OF EMPHASIS: The net trend predicted by the model depends heavily on the specification
of orders within the ARIMA model. It is possible to have several good candidate models with
different predictions of trend. Due to the complexities of ARIMA modeling, the degree of
knowledge needed to interpret the models, and the variability in trends predicted by models all
identified as good candidates, this SOP does not recommend the routine use of ARIMA for
estimating trends.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 21

D� Untitled 8- 1MP
4
0 I
File Edit Tables Rows Cols
Analyze
Graph Tools
View
'I
E1 e 0 � I .. 4-a
A Q
b
giWinndo�wg ®Help
p�
Untitled 8
0
llcdral NOx
NDx Predictian Fermda
Date
Predicted NI
1
1.3
2.9
•
1.3068911663
0112012004
02/10/2004
1.30689116E
Columns (8/0)
2
Actual NOx
A
3
1.1
2.1564910086
03/16/2004
2.15649101
NOx Prediction Formula'
4
1.5
1.7464199244
04/27/2004
1.74641992
Date
5
1.8
1.6692574132
05/25/2004
1.66925741
Predicted NQx
-
6
2.4
1.7175929369
06/15/2004
1.7175929=
.A Std Err Pred NOx
T
1.5
1.93299B8261
07I13#2004
1.93299881
Residual NOx
8
2.7
1.8099906312
08/17/2004
1.8099906=
Upper CL (0.95) NOx
9
1.9
2.0814505252
09/14/2004
2.08145052
Rows
10
0.1
2.03464826
10112121004
2.03464E
All rows
152
11
1.7
1.4703520689
1110912004
1.47035206
Selected
0
12
1.4
1.544%80972
12/07/2004
1.54496809
Excluded
0
13
1.7
1.53968256
01110121005
1.53%82
Hidden
0
141
0.9
1.485246106
0212212005
1.4852461
Labelled
0
14
n 7
1 dd]SR7fi7fi4
III
ngm217fv)r%
1 dd75R7fi7 '
©T
Figure 5. Example Model Prediction Window in JMP
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 22

5.0 Selected Statistical Packages and Statistical Methods
5.1 USGS Kendall Program
The USGS Kendall Program was developed to address a gap in publically available software for
estimating trends using the Seasonal Kendall test and other Kendall tests (Helsel et al., 2006). In
the 1980s USGS popularized Kendall methods, and USGS published computer code supporting
its use in popular statistical packages. However, in later years as statistical analysis moved to
desktop computing, the code became difficult to execute without purchase of commercial
statistical software. As a result, USGS repackaged the code into an executable program which
can be used on computers supporting DOS or DOS emulation. The USGS Kendall program is
freely available and can be downloaded from the following link as of this writing:
http://Pubs.usgs.gov/sir/2005/5275/downloads/
The folder contains the program Kendall.exe, as well as several example input and output files.
Documentation for the program is located here:
http://Pubs.usgs.gov/sir/2005/5275/pdf/sir2005-5275.pdf
The documentation provides an excellent supplemental reference, especially for preparing input
files. The following statistical methods in the USGS Kendall Program are covered by this SOP:
• Mann -Kendall test
• Seasonal Kendall test
• Regional Kendall test
The USGS Kendall Program allows for LOWESS to be used with both the Mann -Kendall and
Seasonal Kendall tests. The documentation does not provide any guidance on selecting f,
although the example output files provided with the program all used a value of 0.5.
Statistics for the Regional Kendall test are identical to those of the Mann -Kendall test. A LOWESS
smooth cannot be used for Regional Kendall, because each site is independent of the other. An
intercept is not reported because, in practice, it would vary by location.
5.2 JMPTM
JMPTM is a commercial statistical software package developed by SAS. It focuses on providing a
powerful yet user-friendly way to view and explore statistical data and conduct analyses. JMP
does not contain the full power or all the analysis options available in SAS, but does cover the
majority of typical user needs. JMP is distributed with varying license levels and add-ons (e.g.,
JMP Pro), but the City's licenses are for the first level, called simply JMP. The SOP covers
statistical methods that can be executed by the first license level. The methods include:
• LOWESS smoothing
• Wilcoxon Rank Sum test
• ARIMA
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 23

While the Wilcoxon Rank Sum test can be performed in JMP, the associated Hodges -Lehman
estimator of step -trend magnitude is not available.
The JMP implementation of ARIMA allows for the full suite of models of the form
ARIMA(p,d,q)x(P,D,Q)S as described in Section 4.5. S is always 12 for monthly data. Models can
be specified with or without a constant (i.e., intercept). JMP does not allow for intervention
analysis within the ARIMA framework, but does provide a full suite of model diagnostic tools
including model fit metrics, ACF, PACF, residuals, and a plot of the predicted model with
confidence limits.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 24

6.0 Procedures
6.1 Preliminary Data Screening
Monitoring data should be screened carefully prior to use in statistical analysis. It is important to
have data with a sufficient period of record for conducting the selected test. Some time series
tests need data to have a consistent sampling frequency, and some tests are more sensitive to
gaps than others. The presence of non -detects complicates some tests, and outliers may also be
problematic. Changing detection limits or severe outliers associated with data entry errors have
the potential to cause spurious identification of trends. The following steps should be performed
for each data set:
When performing trend analyses, knowledge of the system under study is of critical
importance. Random analysis to find trends without a causal explanation risks
identifying a "statistically significant" trend that may not truly exist. It is better to
approach trend analysis with an up -front idea or theory about the outcome. For
instance, it is appropriate to test for trends when management measures have been
implemented in the past upstream of the assessment point. Local knowledge of
watershed conditions should always be used to support findings.
Perform a data count by year to check for the number of annual observations, and to
better visualize gaps and the overall period of record. The YEAR formula in Excel can be
used in a blank column adjacent to the data; next, a Pivot Table can be created that uses
year for summing data counts. Multiple stations or parameters can be queried at one
time in a Pivot Table, making it an efficient tool.
Generally speaking, monitoring data collected on more or less a monthly basis should be
used for trend analysis. It may be appropriate to omit special studies and projects for
two reasons. First, most trend tests require (or are improved by) time series monitored
at a regular interval. Special studies increase the density of data during the period they
are conducted, potentially giving more weight to that time period. Second, samples for
special studies are usually sent to alternate laboratories, which may perform analytical
tests differently (creating bias), or may use different detection limits. Care should also
be taken when mixing data from different projects or laboratories. Staff members
should conduct an extensive comparison of the datasets, checking that the means and
variances are similar. If flow is available, the relationship between concentration and
flow should be compared between the datasets.
• Time spacing between data points should be checked; if there is a relatively high density
of data during some time periods, it may be necessary to select a value representing the
central tendency to allow the data to represent equivalent time periods. The USGS
Kendall Program is less sensitive to this issue because it takes the median of multiple
data points when necessary, but inputs to an ARIMA model should have more or less
equal spacing. If there is more than one observation per month, the data should be
summarized to the monthly level using the median of observations within each month.
• Plot the data versus time to check for obvious outliers; among the tests discussed
ARIMA is the most sensitive to outliers. Also check for changes in detection limits.
Figure 6 shows the pattern for dissolved copper at Third Fork Creek at Hwy 751
(TFO.OTC); the majority of observations were reported at the detection limit until a
change in methods during 2012 resulted in a much lower limit. Figure 7 shows a similar
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 25

problem, where detection limits have declined over time for TP for the Eno River at
West Point (EN8.9ER). Neither data set is very useful for trend analysis because there is
an apparent downward trend associated purely with the change in the detection limit.
Hirsch et al. (1991) note that the presence of less than five percent non -detect records
is not likely to significantly affect the accuracy of a trend slope calculated using one of
the Kendall tests. If the number of non -detects exceeds five percent, the tests can still
be used to determine whether or not there is a trend, but the estimate of the trend
slope becomes unreliable.
• When varying detection limits are present, there are parametric maximum likelihood
estimation methods that can be used to detect trends provided the data are reasonably
normal or can be transformed to normality (Hirsch et al., 1991). Helsel (2010) provides
an excellent review regarding appropriate methods for addressing non -detects,
concluding that maximum likelihood estimation and other methods are superior to
more simplistic assumptions, such as setting the value to one-half the detection limit.
Hirsch et al. (1991) also note that non -parametric methods (i.e., the Kendall tests or the
Wilcoxon Rank -Sum test) can be used if all censored and uncensored values less than or
equal to the highest reporting limit are set equal to the detection limit. However, as
noted in the previous bullet, a large fraction of censored records introduces significant
bias to estimates of trend slope from non -parametric tests.
If the data show a non -normal distribution (almost always the case for fecal coliform
and TSS concentrations, and for loads of any time series), consider whether the data
should be log transformed prior to analysis. For ARIMA, data (after station arization)
should be approximately normal, although there is leeway for some skew to be present.
The Kendall tests and Wilcoxon Rank Sum tests allow for non -normality because the
tests are based on rank, so transformations may not be needed.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 26

8
6
5 • ••• NN••
J
3
4
3
2
Dissolved Copper at Third Fork Creek at Hwy 751
• • NNN NNl•N N
•
•
NN• N•NN••NN•N •
1
01/01/2004 01/01/2006 01/01/2008 01/01/2010 01/01/2012 01/01/2014
Date
Figure 6. Time Series of Dissolved Copper at City of Durham Station TFO.OTC (non -detects
shown at detection limit)
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 27

0.14 TP, Eno river at West Point on the Eno
0.12 •
0.10 • •
0.08 • • • •
E • • • • •
0.06 • N • •« •• •• • • «
•
M«• 69406040D •• •« • • • •
0.04 • •
0.02 N .• • • N
01/01/2004 01/01/2006 01/01/2008 01/01/2010 01/01/2012
Date
Figure 7. Time Series of TP at City of Durham Station EN8.9ER (non -detects shown at detection limit)
6.2 Selection of Appropriate Statistical Tests
After screening the data, the appropriate statistical test(s) is(are) selected based on the
following considerations:
1. If a time series from a single station is being tested for trend, conduct either the Mann -
Kendall test for data that do not show seasonality, or the Seasonal Kendall test for data
that show seasonal variation.
a. See discussion in Section 6.4.2 for guidance on determining whether data show
seasonal variation. Best judgment may be required for borderline cases; always
document assumptions and decisions whenever they are made.
b. If flow data are available at a nearby gaging station, conduct additional Mann -
Kendall or Seasonal Kendall analysis using the LOWESS option in the USGS
Kendall program. Use a range off parameters as discussed in Section 4.1.
Results with LOWESS are preferable to those without because it is likely that
flow explains some of the variation found in concentration data. However, if
load is being analyzed, do not use LOWESS because the data have already been
transformed with flow.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 28

2. If time series from multiple stations are being tested to see if there is an overall trend
across all the stations, use the Regional Kendall test. There is no option for defining
seasonality or using LOWESS with Regional Kendall.
3. If data are being examined to determine if there is a distinct difference in the magnitude
of the values between two time periods, use the Wilcoxon Rank Sum test.
4. Due to the overall complexity of the ARIMA procedure and challenges interpreting
results for trend, ARIMA should not be used for trend analysis in the absence of other
tests. However, ARIMA may have utility for performing a detailed analysis of a time
series for interpretation of trends.
6.3 General Instruction for Documenting Assumptions and Decisions
Whenever judgment is used in conducting analyses, all assumptions and decisions should be
documented clearly with the results by the analyst. Of particular note, significance levels for
critical p values should be agreed upon and documented by staff members prior to conducting
tests.
6.4 Mann -Kendall and Seasonal Kendall Tests
The Mann -Kendall and Seasonal Kendall tests are discussed together, because in practice both
tests should be performed on the same data when seasonality is suspected but not apparent.
The option for LOWESS smoothing using flow is also covered. The first topic is preparation of the
input file, which has very specific requirements. Next, the use of the tests is demonstrated with
City monitoring data for a variety of parameters and locations. The analyses are presented by
pollutant/location, rather than by test. In practice, it is best to conduct different tests on the
same data to better understand the roles of flow and seasonality in explaining portions of
variability in the data.
6.4.1 Preparation of Input Files
The USGS Kendall Program is a Windows executable file run within the DOS environment. To
execute Kendall.exe, double-click on it in the folder where it resides. When the program is run, it
asks for the name of the input file and output file. The input file must be located in the same
folder as Kendall.exe, and the output file will be saved in the same folder. Be sure to use unique
names, and be careful about using the same output filename as one you already have because
the program will overwrite without prompting.
The input file has a very specific fixed -length format; deviations from the format will (at best)
result in a program error, or (at worst) give unexpected or erroneous results. The input file is
prepared in a fixed text format. The first line defines which tests are to be performed, while
subsequent lines contain the input data. The structure of the first line is as follows:
A B CCC DDDDDDDDD...
12345678901234567890
where A, B, C, and D are input variables, and the numbers indicate the specific positions of A
through D. The meanings and possible values for A through D are:
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 29

• A: test to be performed, and structure of data. Placed in column 1.
o A = 1 Seasonal Kendall test; data structure is decimal year, data value, optional
flow value. Season is interpolated from the value of C.
o A = 2 Seasonal Kendall test; data structure is year, season, data value, optional
flow value. The SOP recommends using this method for Seasonal Kendall
because the season is specified rather than interpolated.
o A = 3 Regional Kendall test; data structure is year, site ID, data value.
o A = 4 Mann -Kendall test; data structure decimal year, data value.
• B: Switch for LOWESS option. Placed in column 3.
o B = 0 No LOWESS (test computed on values only).
o B = 1 Use LOWESS (test computed on LOWESS residuals).
• C: Number of seasons, maximum of 120 (only used for A = 1). Placed in column 5-7. Use
a value of 12 to approximate months.
• D: Title printed in output data file (maximum 59 characters). Placed beginning at
column 9.
The data are entered, line by line, below the first line. The following rules should be followed
closely:
• Values can be separated by commas or spaces only — no tabs! This SOP recommends
using commas. Data can be organized in Excel on a separate worksheet and exported as
a CSV (comma separated values) file.
• No blank lines.
• No missing values. If there is a missing value in the input series, delete the line.
• Decimal years have the form 2005.8612, where 2005 is the year and 0.8612 is the
fractional part of the year (corresponding to November 9 in this case). The following
Excel formula can be used to construct the decimal year, where Al contains the date:
=YEAR(Al)+YEARFRAC(DATE(YEAR(Al),l,l),Al,l).
• For Mann -Kendall and Seasonal Kendall tests, flow can be included in the input file even
if it is not used (B = 0). Flow is entered as a second comma -separated value after the
observation on each data line.
• For the Seasonal Kendall test, multiple values within a given year and season are
allowed. The USGS Kendall program calculates the median of the values in keeping with
the procedure used by the original Seasonal Kendall test.
An example input file is shown here for a Seasonal Kendall Test (A = 2) with LOWESS on flow (B =
1). CCC is not used because A is 2 and not 1. The A B CCC DDDDD... line is shown above the first
line for reference (it is not part of the file). In this case, the data have four values — year, season,
NOx, and flow. The flow units are cfs, although the units are not germane so long as the same
units are used for all measurements.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 30

A B CCC DDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD
2 1 Seasonal Kendall with LOWESS, NOx at EN8.9ER
2004,1,0.8,57
2004,2,0.6,130
2004,3,0.5,186
2004,4,0.6,59
2004,5,0.5,21
2004,6,0.4,17
2004,7,0.1,2.1
2004,8,0.4,278
2004,9,0.6,55
2004,10,0.3,133
2004,11,0.1,36
2004,12,0.6,71
2005,1,0.5,55
2005,2,0.5,69
2005,3,0.5,396
2005,4,0.2,146
2005, 5, 0 .5, 37
2005, 6, 0.3, 49
2005,7,0.2,7.2
6.4.2 Example: TP at City of Durham Station EL1.9EC
Trends in TP at EL1.9EC (Ellerbe Creek at Glenn Rd.) are investigated using Mann -Kendall and
Seasonal Kendall tests. LOWESS smoothing using flow is presented as well. As discussed
previously in Section 6.1, screening should be conducted to ensure the data are suitable for time
series analysis. In the case of EL1.9EC, ambient monitoring data have been collected
continuously since 2004, while flow monitoring at the co -located USGS station 02086849 began
in 2006. TP monitoring shows no significant gaps. A plot of TP versus time (Figure 8) shows no
major reductions in in the detection limit over time. Approximately ten percent of the data are
reported as being below a detection limit. As discussed in Section 6.1, as a general rule there
should be no more than five percent of the samples reported below the detection limit or the
Kendall tests risk bias in the estimate of trend slope. However, the majority of observations are
considerably higher than the detection limit, so uncertainty in the lowest values is less likely to
be an issue. The values were set equal to the reported detection limits in subsequent analyses.
Note that data from any special studies were removed from the analysis to minimize issues with
the use of different analytical laboratories and/or methods — only ambient monitoring data
were used. TP data collected prior to the inception of flow monitoring in 2006 were removed
from the analysis so that tests with and without use of LOWESS for correcting for correlation to
flow would be comparable. Because a normal distribution assumption is not a requirement for
the Kendall tests, the data were not log transformed. First, the Mann -Kendall tests are
performed. Next, the issue of seasonality of TP is explored, and the tests are repeated using the
Seasonal Kendall test.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 31

TP at EL1.9EC, Ellerbe Creek at Glenn Rd.
1.4
•
1.2
1.0
J
0.8
•
LCJ1
E
d
H
0.6
•
0.4
•
•
0.2
•
b • '
M• • •
0.0
01/01/2004
01/01/2006 01/01/2008 01/01/2010 01/01/2012 01/01/201
Date
Figure 8. Time Series of TP at City of Durham Station EL1.9EC
First, an input file is prepared for the Mann -Kendall test. The first few lines are shown:
4 0 MK test, TP at EL1.9EC
2006.024658,0.14,13
2006.10137,0.08,18
2006.19726,0.09,17
Note that the switch for LOWESS is set to 0 (no LOWESS) for the first test. Flow values are
present in the input file, but are simply ignored. The Kendall.exe program is executed in a folder
including the input file. The program asks for the filename of the input file, then asks for an
output filename. Even though Kendall.exe is running under DOS emulation, input and output
filenames can be longer than eight characters. The following output is produced:
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 32

Kendall's tau Correlation Test
US Geological Survey, 2005
Data set: MK test, TP at EL1.9EC
The tau correlation coefficient is 0.236
S = 1008.
z = 3.357
p = 0.0008
The relation may be described by the equation:
Y =-14.888 + 0.7452E-02 * X
The title included in the input file is reported in the output file. The reported statistics were
previously discussed in Section 4.2. Note that 5 is a large positive number indicating a positive
trend, with a corresponding p showing a high level of significance. The magnitude of the trend
found using Sen's non -parametric slope estimator (Sen, 1968) is reported as approximately-
0.0075 mg/L per year. The mean of TP during the monitored period is 0.125 mg/L so the trend is
weak but significant according to the Mann -Kendall test.
Next, the test is run with LOWESS to account for possible impacts of correlation between
concentration and flow. The first few lines of the input file are as follows:
4 1 MK test with LOWESS (f = 0.25), TP at EL1.9EC
2004.054645,0.8,57
2004.112022,0.6,130
2004.20765,0.5,186
Note the only difference is that the switch for LOWESS is set to 1, and the title is altered. Flow
values are now processed. The Kendall.exe program is executed, and input/output filenames
specified. Next, the program asks for the value of the f parameter for LOWESS. Document f in
the input file title as shown above to assist with tracking output results. Run at least three
iterations of the Kendall program when using LOWESS, with f set to 0.2, 0.5, and 0.8, based on
the range of typical f values recommended by Cleveland (1979). Adjust f as needed if the Kendall
program reports that it changed f to 1.0, which may occur if f is too small for the dataset (in this
example, f = 0.25 was used as a value of 0.2 was too small for the dataset). An examination of
the results should provide a good indication whether the magnitude off influences the
outcome. The results for f = 0.5 are shown:
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 33

Kendall's tau Correlation Test
US Geological Survey, 2005
Data set: MK test with LOWESS (f = 0.5), TP at EL1.9EC
Flow -adjusted concentrations were calculated as
residuals from a LOWESS smooth line (f = 0.50)
for concentration vs. streamflow
The tau correlation coefficient is 0.222
S = 951.
z = 3.153
p = 0.0016
The relation may be described by the equation:
Residual =-14.388 + 0.7155E-02 * X
The value of p remains highly significant, and the reported trend slope is very similar to the test
without LOWESS, indicating that the trend observed without accounting for flow was a valid
trend. For f = 0.25 the value of p is 0.0022 and the trend is 0.0065, while for f = 0.8, the value of
p is 0.0022 and the trend is-0.0073. Even with a wide range off the results converge on a
positive trend in the neighborhood of 0.007 mg/L per year.
It is possible that TP shows seasonal variation. If so, it is difficult to discern using the time series
plot in Figure 8. The user should plot the variations of the observations relative to the mean by
month to see if seasonality is evident; Figure 9 includes the median among all the months as a
red square to indicate central tendency. There appears to be a decrease in TP during spring
months relative to summer months, although there is a fair amount of scatter.
TP at EL1.9EC, 2006 - 2013
Monthly Variation Against Mean (red = monthly median)
1
0.8
0.6
0.4
0.2
0
-0.2 1 1
1 2 3 4 5 6 7 8 9 10 11 12
Figure 9. Monthly Variation of TP at City of Durham Station EL1.9EC from 2006 - 2013
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 34

Another useful way to diagnose seasonality is by using the ACF introduced in Section 4.5. This
can be obtained in JMP by running Analyze > Modeling > Time Series, specifying TP as the input
for Y, Time Series, and date as X, Time ID. The resulting ACF is shown in Figure 10 on the left.
Seasonality is identified by pattern of alternating positive and negative lags at distances of 12
periods. The signature of seasonality in the ACF is very weak and not statistically significant,
consistent with the monthly variation plot in Figure 9.
Time Series TP
1
Mean
0.1260778
0.8
Std
0.1395479
N
90
a 0.6
Zero Mean ADF
-5.261318
~ 0.4
Single Mean ADF
-8.087721
Trend ADF
-8.672671
0.2
0
01/01/2004 01/01/2007
01/01/2010 01/01/2013
Date
Time Series Basic Diagnostics
Lag
AutoCorr-.8-.6-.4-.2 0 .2 .4 .6 .8 Ljung-Box Q
p-Value
Lag Partial-.8-.6-.4-.2 0 .2 .4 .6 .8
0
1.0000
0 1.0000
1
0.1222
1.3896
0.2385
1 0.1222
2
0.1253I
2.8675
0.2384
2 0.1121
3
0.0083
2.8740
0.4115
3-0.0196
4
-0.0327
2.9770
0.5617
4-0.0469
5
-0.0468
3.1907
0.6706
5-0.0378
6
0.0031
3.1917
0.7844
6 0.0228
7
-0.0023
3.1922
0.8667
7 0.0055
8
-0.0170
3.2215
0.9197
8-0.0229
9
-0.0849
3.9581
0.9141
9-0.0875
10
-0.0162
3.9855
0.9480
10 0.0062
11
0.0325
4.0964
0.9670
11 0.0590
12
0.1189
5.5967
0.9350
12 0.1149
13
0.1399
7.7026
0.8624
13 0.1013
14
-0.0121
7.7184
0.9035
14-0.0791
15
0.1758
11.1314
0.7432
15 0.1682
16
-0.0168
11.1630
0.7993
16-0.0299 . .
.
17
0.0032
11.1641
0.8479
17-0.0156
18
0.0000
11.1641
0.8873
18 0.0010
19
0.2399
17.8756
0.5308
19 0.2666
20
0.0161
17.9063
0.5936
20-0.0201
21
0.0650
18.4132
0.6227
21 0.0262
22
0.1186
20.1271
0.5751
22 0.1461
23
0.0604
20.5780
0.6068
23 0.0481
24
0.0659
21.1224
0.6315
24 0.0663
25
-0.0613
21.6014
0.6586
25-0.1306
Figure 10. ACF for TP at City of Durham Station EL1.9EC
Based on these analyses, it unlikely that seasonality of statistical significance is present in the TP
data. However, the Seasonal Kendall test is performed for demonstration purposes. An input file
is prepared for the Seasonal Kendall test, first without LOWESS:
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 35

2 0 SK test,
2006,1,0.14,13
2006,2,0.08,18
2006,3,0.09,17
TP at EL1.9EC
Note that dates are now simply reported using year, but season is included after date on the
input rows. In this example, seasons are defined using month. Kendall.exe is executed, and
input/output filenames specified, yielding the following results:
Seasonal Kendall Test for Trend
US Geological Survey, 2005
Data set: SK test, TP at EL1.9EC
The record is 9 complete water years with 12 seasons per year
beginning in water year 2006.
The tau correlation coefficient is 0.283
S = 89.
z = 3.336
p = 0.0008
p = 0.0251 adjusted for correlation among seasons
(such as serial dependence)
The adjusted p-value should be used only for data with
more than 10 annual values per season.
The estimated trend may be described by the equation:
Y = 0.4500E-01 + 0.1000E-01 * Time
where Time = Year (as a decimal) - 2005.75 (beginning of
first water year)
Two p values are reported; the second includes the correction for serial correlation among
values across years. The data do not quite meet the requirement for using the adjusted p value
(at least 10 years of data distributed across 12 months). One aspect of the results is important to
note — when serial correlation is considered, p increases to a much higher value (indicating
weaker statistical significance). Additional runs were conducted with LOWESS using the series of
three f values. The p values and trends are shown in Table 2, along with the results of all the
previous tests for reference in the discussion that follows.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 36

Table 2. All Kendall Test Results for TP at City of Durham Station EL1.2EC
Analysis
p
Adjusted p
Trend
(mg/L/yr)
Mann -Kendall, no LOWESS
0.0008
0.0075
Mann -Kendall with LOWESS, f = 0.25
0.0022
0.0065
Mann -Kendall with LOWESS, f = 0.5
0.0016
0.0072
Mann -Kendall with LOWESS, f = 0.8
0.0022
0.0073
Seasonal Kendall, no LOWESS
0.0008
0.0251
0.0100
Seasonal Kendall with LOWESS, f = 0.25
0.0007
0.0166
0.0083
Seasonal Kendall with LOWESS, f = 0.5
0.0003
0.0119
0.0097
Seasonal Kendall with LOWESS, f = 0.8
0.0002
0.0088
0.0107
At first glance these results seem difficult to interpret. There is no clear indicator of which test
best represents the process. Are the results with LOWESS on flow preferable? If so, which f
value is appropriate? Should the Seasonal Kendall test results be used instead of the Mann -
Kendall test results? The answer comes down to best judgment, but the underlying question
should also be considered — Is there an upward trend in TP, and if so what is its significance and
magnitude? First, the Mann -Kendall results without and with LOWESS can be examined. Adding
LOWESS does not improve the significance of p, suggesting that flow does not explain a
significant portion of the variance. Trends estimated from the different analyses are very close
to each other, so there is consistency among the different methods. Second is the question of
Seasonal Kendall without and with LOWESS. As shown previously, the data do not appear to
show seasonality, but with the addition of LOWESS, the Seasonal Kendall tests perform
somewhat better than the Mann -Kendall tests based on the unadjusted p values, although only
marginally. However, the adjusted p values for Seasonal Kendall are higher than the Mann -
Kendall values, and with eight years of data the time series comes close to the ten year
threshold recommended for using the adjusted p values. The range in Seasonal Kendall trends is
more divergent than for the Mann -Kendall tests, but all show higher trend with time than the
Mann -Kendall results. Given the lack of significant seasonality in the data, the Mann -Kendall
results should be favored, and the trend in TP is in the neighborhood of 0.007 mg/L per year. As
noted in the introduction, whenever judgment is used in conducting analyses, all assumptions
and decisions should be clearly documented.
6.4.3 Example: Water Quality Index at EN8.9ER
The Water Quality Index (WQI) is a relative measure of water quality reported by the City. It is
calculated as a composite from several parameters, and is scaled from zero up to 100 (signifying
excellent quality). Ambient monitoring WQI values at EN8.9ER (Eno River at West Point) have
been calculated continuously beginning in 2004. Monitoring shows no significant gaps, although
values were not calculated for seven of the months during the ten year period. A plot of WQI
versus time (Figure 11) shows high values for the majority of observations. The data visually
suggest a gradual increase in WQI, although the trend may not be significant given the variability
in the observations. Note that data from any special studies were removed from the analysis to
minimize issues with the use of different analytical laboratories and/or methods — only ambient
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 37

monitoring was used. It is possible that WQI has a seasonal trend. A plot of variations of the
observations relative to the mean by month is shown in Figure 12; the medians of the variations
are skewed towards high values. There does not appear to be any seasonal pattern. This
conclusion is reinforced by the ACF (Figure 13); not only is there no indication of seasonal lag,
but the entire series has no significant lags or patterns in the lags, and appears to be stationary
(supported by the series of Ljung-Box Q p-values that show no significance). Given the lack of
seasonality, Seasonal Kendall tests are not conducted; only Mann -Kendall tests are performed.
Figure 11. Time Series of WQI at City of Durham Station EN8.9ER
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 38

10
5
0
-5
-10
-15
-20
-25
-30
-35
-40
WQI at EN8.9ER, 2004 - 2013
Monthly Variation Against Mean (red = monthly median)
1 2 3 4 5 6 7 8 9 10 11 12
Month
Figure 12. Monthly Variation of WQI at City of Durham Station EN8.9ER from 2004 — 2013
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 39

Time Series WQI
Mean
93.554545
100
"
•ifs
Std
8.84779
90
N
110
r
Zero Mean ADF-0.407759
80
Single Mean
ADF-10.95338
70
Trend ADF
-11.20664
60
01/01/2004 01/01/2007
�r-
01/01/2010 01/01/2013
Date
Time Series Basic Diagnostics
Lag
AutoCorr-.8-.6-.4-.2 0 .2 .4
.6 .8 Ljung-Box Q
p-Value
Lag
Partial
8-.6-.4-.2 0 .2 .4 .6 .
0
1.0000
0
1.0000
. . '
1
-0.02851
0.0918
0.7618
1
-0.0285
2
-0.0319
0.2076
0.9014
2
-0.0327
3
-0.0056
0.2112
0.9758
3
-0.0075
4
-0.0943
1.2444
0.8707
4
-0.0959
5
0.1821
:7 5.1343
0.3997
5
0.1779
6
-0.0050
5.1372
0.5263
6
-0.0038
7
0.0957
6.2328
0.5128
7
0.1113
8
-0.0362
6.3911
0.6035
8
-0.0432
9
-0.1092
7.8444
0.5499
9
-0.0716
10
-0.0148
7.8713
0.6414
10
-0.0579
11
0.0757
8.5848
0.6602
11
0.0930
12
-0.0003
8.5848
0.7379
12
-0.0474
13
-0.1299
10.7291
0.6335
13
-0.1315
14
0.0191
10.7761
0.7035
14
0.0303
15
-0.0286
10.8823
0.7609
15
-0.0030
16
-0.0676
11.4815
0.7788
16
-0.0886
17
0.0100
11.4947
0.8297
17
-0.0077
18
-0.08991
12.5762
0.8161
18
-0.0775
19
-0.1469 E
15.4978
0.6905
19
-0.1773
20
0.0244
15.5793
0.7424
20
-0.0029
21
0.0033
15.5808
0.7927
21
0.0080
22
-0.0245
15.6650
0.8322
22
-0.0849
23
-0.0820
16.6182
0.8276
23
-0.0820
24
0.0198
16.6742
0.8622
24
0.0953
25
0.0819L17.6456
0.8571
25
0.0755
Figure 13. ACF for WQI at City of Durham Station EN8.9ER
Input files were prepared for the Mann -Kendall test without LOWESS, and with LOWESS using
f = 0.2, f = 0.5, f = 0.6, and f = 0.8. See Section 6.4.1 for instructions regarding preparation of
input files, and Section 6.4.2 for example input file excerpts. Results are shown in Table 3. S
values are positive, indicating an upward trend - in other words, an improvement in the WQI.
Addition of LOWESS improves p across a range off, so it appears that flow explains part of the
variation in WQI. p values are similar across the range off. What's interesting is the prediction of
trend does not change linearly with variation in f. In fact, there is a drop in the estimated trend
from 0.2414 for f = 0.5 to 0.2061 for f = 0.6. Given the variation in trend with changes in f, the
trend is reported as a range rather than attempting to select one value - i.e., there is an
increase in WQI of about 0.20 - 0.25 per year at EN8.9ER with statistical significance.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 40

Table 3. Mann -Kendall Tests for NOx at City of Durham Station EN8.9ER
Analysis
S
p
Trend
(change/yr)
Mann -Kendall, no LOWESS
760
0.0473
0.1487
Mann -Kendall with LOWESS, f = 0.2
840
0.0302
0.2551
Mann -Kendall with LOWESS, f = 0.5
838
0.0306
0.2414
Mann -Kendall with LOWESS, f = 0.6
838
0.0306
0.2061
Mann -Kendall with LOWESS, f = 0.8
818
0.0348
0.2072
6.4.4 Example: Dissolved Oxygen at City of Durham Station NH3.3SC
Ambient monitoring of DO concentrations at NH3.3SC (Sandy Creek at Cornwallis Rd.) has been
conducted continuously beginning in 2008. Monitoring shows no significant gaps, and there
were no values below detection limits. A plot of DO versus time shows a strong seasonal pattern
in measurements (Figure 14), which is reinforced by the strong seasonal pattern visible in the
ACF (Figure 15). Note that data from any special studies were removed from the analysis to
minimize issues with the use of different analytical laboratories and/or methods. Monitoring
from 2011 on suggest a possible upward trend. Mann -Kendall tests are not appropriate given
the strong seasonal signature, so only Seasonal Kendall tests are performed. Flow monitoring
began at the site in August of 2008, so tests with LOWESS are also executed, noting that the
period of record with LOWESS is truncated seven months relative to the test without LOWESS.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 41

Figure 14. Time Series of DO at City of Durham Station NH3.3SC
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 42

Time Series DO
14
Mean
8.0495714
12
Std
2.3723614
N
70
10
O
Zero
Mean ADF
-1.224824
8
Single Mean ADF
-3.876239
6
Trend ADF
-3.857741
4
01/01/2007 01/01/2009 01/01/2011
01/01/2013
01/01/2015
Date
Time Series Basic Diagnostics
Lag
AutoCorr-.8-.6-.4-.2 0 .2 .4 .6 .8
Ljung-Box Q
p-Value
Lag
Partial-.8-.6-.4-.2
0 .2.4.6.
0
1.0000 :
0
1.0000 1
1
0.6276
28.7672
<.0001*
1
0.6276
2
0.3524
37.9694
<.0001*
2
-0.0684
3
-0.0249
38.0160
<.0001*
3
-0.3617
4
-0.3728 1
48.6302
<.0001*
4
-0.3644
5
-0.5953
76.1130
<.0001*
5
-0.2740
6
-0.7124
116.085
<.0001*
6
-0.3354
7
-0.5556
140.780
<.0001*
7
-0.0702
8
-0.3615
151.404
<.0001*
8
-0.2031
9
-0.0273
151.465
<.0001*
9
-0.0682
10
0.2693
157.556
<.0001*
10
-0.0654
11
0.5408
182.542
<.0001*
11
0.1142
12
0.6650
220.971
<.0001*
12
0.1303
13
0.6031
253.137
<.0001*
13
0.1343
14
0.3222
262.481
<.0001*
14
-0.1583
15
0.0212
262.522
<.0001*
15
-0.1331
16
-0.2784
269.757
<.0001*
16
0.0624
17
-0.5040
293.912
<.0001*
17
0.0436
18
-0.5692
325.314
<.0001*
18
0.0059
19
-0.4756
347.669
<.0001*
19
0.0602
20
-0.2969
356.555
<.0001*
20
-0.0916
21
-0.0436
356.750
<.0001*
21
-0.0480
22
0.1994
360.925
<.0001*
22
-0.0329
23
0.3830
376.652
<.0001*
23
-0.1000
24
0.4636
400.198
<.0001*
24
-0.1197
25
0.4571
423.595
<.0001*
25
0.0297
Figure 15. ACF for DO at City of Durham Station NH3.3SC
Input files were prepared for the Seasonal Kendall test without LOWESS, and with LOWESS using
f = 0.2, f = 0.5, and f = 0.8. See Section 6.4.1 for instructions regarding preparation of input files,
and Section 6.4.2 for example input file excerpts. Results are shown in Table 4. This time the S
values are small and the corresponding p values large (all much greater than 0.05), indicating
that the null hypothesis of no trend cannot be rejected. One interesting question relates to how
to interpret the default p versus p adjusted for autocorrelation. The data count does not meet
the recommended minimum of ten years of monthly data to use the adjusted p, which suggests
the default p is appropriate. Hirsch et al. (1982) recommend using a minimum of three years of
monthly data to conduct the Seasonal Kendall test; in this example, there are about five years,
so the default p is likely appropriate. The reason for the ten year requirement for the adjusted p
given by Hirsch and Slack (1984) is that the calculation of the covariance matrix to arrive at the
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 43

adjusted p becomes inexact for fewer than about 120 observations. But what if there are nine
years of data, and the default p and adjusted p are markedly different from each other? Best
judgment should be used, but it seems reasonable that as the number of years of monthly
monitoring approaches ten, the true p likely becomes closer to the adjusted p. Regardless,
assumptions and decisions should always be documented with the results.
Table 4. Seasonal Kendall Tests for DO at City of Durham Station NH3.3SC
Analysis
S
p
Adjusted
p
Trend
(change/yr)
Seasonal Kendall, no LOWESS
14
0.4651
0.4294
0.0550
Seasonal Kendall with LOWESS, f = 0.2
5
0.7941
0.8009
0.0869
Seasonal Kendall with LOWESS, f = 0.5
15
0.3611
0.4228
0.1391
Seasonal Kendall with LOWESS, f = 0.8
19
0.2403
0.4032
0.1813
6.5 Regional Kendall Test
General guidance on preparation of input files for using the USGS Kendall program is provided in
Section 6.4.1.
The City has implemented a number of management measures in the Ellerbe Creek watershed
in recent years. Have conditions improved overall in Ellerbe Creek? The Regional Kendall test
can be used to examine the possibility of an overall trend across multiple stations. To
demonstrate the use of the Regional Kendall test, fecal coliform bacteria from multiple Ellerbe
Creek watershed stations are analyzed. Station locations are shown in Figure 16. First, the data
are screened to look for overlapping periods of record (Table 5). Note that data from special
studies were removed from the analysis to minimize issues with the use of different analytical
laboratories and/or methods — only ambient monitoring was used.
When selecting stations for the analysis, the goal is to maximize both the number of stations
tested and the number of years included to gain the most statistical power, and to avoid
including locations that provide largely redundant information because they are within a short
distance on the same stream. First, the station locations are screened to check for those with
close proximity to each other. The following cases are noted:
• EL5.OEC and EL5.6EC on Ellerbe Creek are within 0.6 miles of each other.
• EL7.9EC and EL8.2EC on Ellerbe Creek are within 0.3 miles of each other.
• EL7.ISEC is downstream of the confluence of South Ellerbe Creek and a tributary of the
creek, and may reflect a composite of conditions at EL7.6SECT and EL8.5SEC.
In the first two cases, the choice of stations to eliminate is easy, because monitoring ended at
EL5.OEC in 2005, and monitoring at EL8.2EC is sporadic. For the third case, best judgment must
be used. One could argue that because EL7.ISEC is about 1.4 miles downstream of EL8.5SEC that
it is different enough to be retained. However, to be conservative EL7.ISEC was removed from
the analysis, a decision helped in part by lack of monitoring in 2004 and 2012. Next, other
stations with significant monitoring gaps are removed. Several of the stations were not
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 44

monitored in 2012, but recall that gaps are allowed in the Kendall tests, so these are retained.
Monitoring began at EL8.6SECUT in 2008, too short a time period to be realistically paired with
the other stations. Eight stations remain for analysis. Among all the data, less than five percent
of the observations were reported as above or below reporting limits. As always, any decisions
and judgment calls should always be documented by the analyst.
Legend
Durham County r
Durham City �01101
A Ellerbe Creek Stations
r
l
ELIOJEC I ss e
nt
Club j +
o,
0
.2EC
EL7.9EC��
EL5;OEC1,
i L,EI
j EL5.6EC
IC E'LU7aeo
~ Gam%
'EL7.6SEICT"
go
if
�.�. •�
`� I �� s
m
m
ham c--
,
f
�
ary
nity
�IJRHAAA
�3��
i
1
r9ryQEL8:6SECUT
d O
r�
G
cna
J� r
I u a� faEL8.1
GC
Q /
oreh a
�¢
/ b �
0 0.5 1
1869
Fast �`�
�`ty O
mommom== Miles
y
Durham Trend Analysis SOP
❑
TETRA
TECH
Ellerbe Creek Monitoring Stations
Information depicted hereon is for reference purposes only and is compiled from the best available sources.
The City of Durham/Durham County assumes no responsibility for errors arising from use or misuse of this map.
Map produced by S. Job, 01-22-2014
NAD 1983 State Plane North Carolina FIPS 3200 Feet
Figure 16. Durham Monitoring Stations in the Ellerbe Creek Watershed
Table 5. Fecal Coliform Data Count for Ellerbe Creek Stations
Station
2004
2005
2006
2007
2008
2009
2010
2011
2012
2013
EL1.9EC
12
12
12
10
12
12
11
12
12
12
EL10.7EC
12
12
12
10
12
12
11
12
12
12
EL5.OEC
12
1
EL5.5GC
12
12
12
9
12
12
11
12
12
EL5.6EC
12
12
12
10
12
12
11
12
12
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 45

Station
2004
2005
2006
2007
2008
2009
2010
2011
2012
2013
EL7.1SEC
11
12
10
12
12
11
12
12
EL7.6SECT
12
12
12
10
12
12
12
12
12
EL7.9EC
12
12
12
10
12
12
11
12
12
12
EL8.1GC
12
12
12
6
12
12
11
12
12
EL8.2EC
11
12
11
10
1
12
12
EL8.5SEC
12
12
12
10
12
12
11
12
12
EL8.6SECUT
12
12
12
12
12
12
An input file is prepared for the Regional Kendall test for the eight locations for a time period
spanning 2004 through 2013. The format is similar to that of the Seasonal Kendall test — only
year is represented — but instead of season in the second data column, a location number is
entered. As noted previously, there is no option for LOWESS in the Regional Kendall test.
Because a normal distribution assumption is not a requirement for the Kendall tests, the data
were not log transformed. The first few lines are shown. Recall the first input variable specifies
the test to be run — in this case, 3 is used for the Regional Kendall test:
3 Ellerbe Creek RK, fecal coliform
2004,1,81
2004, 1, 62
2004,1,6000
2004,1,6000
Note that the data snippet shown above represent values from January through April 2004 at
the station with id = 1. Data sequence within a year is disregarded in Regional Kendall; the test
looks at year to year changes across the population of stations, and detects the overall annual
trend. The following output is produced:
Regional Kendall Test for Trend
US Geological Survey, 2005
Data set: Ellerbe Creek RK, fecal coliform
The record is 10 years at 8 locations
beginning in year 2004.
The tau correlation coefficient is -0.473
S = -149.
z = -5.122
p = 0.0000
The estimated median trend throughout the region
during years 2004 through 2013 is:
Change in Y = -35.42 per year.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 46

The results indicate an overall trend across all the stations with a very high level of significance
(p < 0.0001). S is negative indicating a downward trend, which is reported to be a reduction of
about 35 cfu/100ml- per year.
It is worth considering the impact of an entire year of data (2012) missing from five of the eight
stations. An alternate run is performed with monitoring from 2004 through 2011. Drawbacks
include loss of statistical power due to having eight years of data rather than ten, as well as
dropping the two most recent years of monitoring. The results of the run with reduced period of
record follow:
Regional Kendall Test for Trend
US Geological Survey, 2005
Data set: Ellerbe Creek RK, fecal coliform (short POR)
The record is 8 years at 8 locations
beginning in year 2004.
The tau correlation
S = -70.
z = -3.018
p = 0.0025
coefficient is -0.312
The estimated median trend throughout the region
during years 2004 through 2011 is:
Change in Y = -32.28 per year.
The outcome is essentially the same. p is now measureable but still very small, and the
magnitude of the trend is within ten percent of the run with the longer period of record.
So what does it mean to have an annual overall reduction of fecal coliform of 35 cfu/100 mL per
year? If the typical concentration is 100 cfu/100 mL, then an annual reduction of 35 would
indeed be remarkable. However, if the typical concentration is 2,000 cfu/100 mL, then a
decrease of 35 is less than a two percent annual decline. Recall that the trend reported by the
Kendall tests, Sen's non -parametric slope estimator (Sen, 1968), represents the median value
from the population of pairwise trends, so it stands to reason that the central tendency of the
pooled dataset is a good indication of a benchmark to which the reduction can be compared.
Because fecal coliform data are typically non -normal (definitely the case for Ellerbe Creek
measurements), then a non -parametric indicator of central tendency is appropriate. The
geometric mean and median of the pooled dataset are respectively 480 and 430 cfu/100mL, so a
decrease of 35 cfu/100mL represents an annual reduction in the neighborhood of 8 percent.
6.6 LOWESS
JMP provides a way to perform LOWESS, although it is referred to by another name ("Local
Smoother") within the JMP menu structure. In addition, one can save the residuals from a
LOWESS smooth to perform subsequent analyses. To execute LOWESS, perform the following
steps:
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 47

1. Ina data window, select Analyze > Fit Y by X.
2. Select X, Factor, and Y, Response variables (for example, flow for X and turbidity for Y).
Click OK. The Bivariate Fit window appears (Figure 17).
3. Click the red triangle pointing down next to the text "Bivariate Fit of yyyy by xxxx"
(substitute variable names for yyyy and xxxx). Select Kernel Smoother.
4. Change Smoothness (alpha) to the desired value. Note that this parameter is equivalent
to f in the USGS Kendall Program.
5. Change Robustness to the desired number of passes. The SOP recommends no more
than 2 passes, based on the recommendation of Cleveland (1979).
6. To save the residuals, select the red triangle pointing down next to "Local Smoother",
and click Save Residuals. The residuals are now available in the original data window.
As an example, the residuals from a LOWESS smooth of turbidity at TFO.OTC using flow is shown
in Figure 18. The reduction in variability is apparent when compared to the original turbidity
series, shown as individual points in the figure. Note that flow monitoring began in August 2008,
so values from earlier in 2008 are truncated from the LOWESS smoothing operation and
resulting turbidity residuals.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 48

Bivariate Fit of Turb. By Flow
250
200
150
100
50
•
•
% •
0
0.001 0.01
0.03 0 1 0.2 0.4 1 2 3 4 10 20 40 100 300 600
Flow
—Local Smoother
Local Smoother
R-Square
0.242791
Sum of Squares Error
85609.86
Local Fit (lambda)
Linear
Weight Function
Tri-cube
Smoothness (alpha):
0.5
Robustness
2 passes
Figure 17. Bivariate Fit Window in JMP, Showing Options for Kernel Smoother
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 49

Turb. & Residuals Turb. at TFO.OTC, Third Fork Creek at Hwy
751
200
150
100
•
• • • •
• •
0
200
150
•
100
50
v
0
•
• •: ..••••• ..••••• •••••••••• •••
•
.•. •••.•• •.
-50
•
-100
1/1/2008
1/1/2009 1/1/2010 1/1/2011 1/1/2012
1/1/2013 1/1/2014
Date
Figure 18. Turbidity and Turbidity Residuals from LOWESS Smooth using Flow, City of Durham
Station TFO.OTC
6.7 Wilcoxon Rank Sum Test
The Wilcoxon Rank Sum test is useful for testing whether two populations have different means.
For example, the Wilcoxon test can be used to investigate whether a step change in observed
concentrations exists at a statistically significant level, or, when there has been a wide gap in
monitoring, one may wish to test whether the early monitoring data are from the same
population as the later monitoring data. Environmental monitoring data frequently do not meet
the requirement of normality for the independent sample t test. The non -parametric Wilcoxon
test is not subject to this limitation.
For example, the Triangle WWTP performed an upgrade in technology between April and May
of 2005, adding a 5-stage biological nutrient removal system. The Upper Cape Fear River Basin
Association (UCFRBA) conducted monitoring downstream of the plant outfall from 2000 through
2012. The data provide a way to test whether pollutant concentrations have changed following
the upgrade. Use of the Wilcoxon Rank Sum test is demonstrated using JMP.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 50

First, the data were preprocessed in Excel. Data sets were prepared for TN, TP, and turbidity.
TKN and NO2+NO3 were summed (when both were reported) to calculate TN. Because the
Wilcoxon Rank Sum test uses rank rather than value, measurements at the detection limit were
set equal to the detection limit as a conservative assumption. Less than five percent of the TN
values were affected by detection limits, and none of the TP values were reported below the
detection limit. Because the Wilcoxon Rank Sum test is based on testing populations from two
different time periods, a categorical variable is needed to conduct the analysis in JMP. An
additional column named "Group" is added to categorize the observations into "Pre" versus
"Post," with the dividing line set between April and May of 2005, corresponding to the timing of
the plant upgrade. Plots of TN, TP, and turbidity (log transformed) are shown in Figure 19. There
is a clear visual reduction in TN, but the possible presence of any step change for TP and
turbidity is difficult to discern.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 51

TN at B3670000, Northeast Creek at O'Kelly Church Rd
30
25
20
J
ego
1.s
1.6
1.4
1.2
J
bD 1
E
o.s
0.6
0.4
0.2
0
M
TP at B3670000, Northeast Creek at O'Kelly Church Rd
M
s♦ ♦
2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012
loot
10(
Turbidity at 63670000, Northeast Creek at O'Kelly Church Rd
0
0
1
2000 2001 2002 2003 2004 2005 2006
2007 2008
2009 2010 2011 2012
Figure 19. Time Series of TN, TP, and Turbidity at UCFRBA Station B3670000,
Downstream of Triangle WWTP
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 52

To perform the Wilcoxon Rank Sum test in JMP, undertake the following steps:
1. Ina data window, select Analyze > Fit Y by X.
2. Select X, Factor, and Y, Response variables (see Figure 20). The X variable must be
categorical. In this example, TN is selected for Y. Click OK.
3. Optionally, click the red triangle pointing down next to the text "Bivariate Fit of yyyy by
xxxx" (in this example, "TN by Group"). Select Display Options > Box Plots.
4. Click the same red triangle again and select Nonparametric > Wilcoxon Test.
5. The output of the test is shown below the plot (Figure 21). Use the 2-Sample Test,
Normal Approximation results to test the null hypothesis Ho that the means of the
populations are the same. In this example, Ho is rejected and we can say that the means
are indeed different, as supported by the high significance of the p value at <0.0001.
Mean TN concentrations before and after the upgrade are 9.05 and 2.19 mg/L,
respectively, while median TN concentrations before and after the upgrade are 7.8 and
1.74 mg/L, respectively.
Fit Y oV X- Contextual - JMP _ y b 4 0 I 0
Distribution ofYfor each X. Modeling types determine analysis.
Select Columns —
5 Columns
................................
AD ate
iTurbirlity
'DTP
:Group
.................................................................................
0newa
A
IL
J
001
Bivariate
Oneway
11
Logistic
Contingency
Ad Cd
Cast Selected Columns into Roles
Y, Responsel JTN
bona t
Action
OK
Cancel
Factar� d.GrvuP Remove
cp bona t
Re€all
Help
Black aptionai
Weight optionm numeric
Freg aptionatnumerir
By optional
Figure 20. Bivariate Window in JMP for Wilcoxon Rank Sum Test
i ❑r'
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 53

Oneway Analysis of TN By Group
25 •
20
z 15
10
5
0
Post Pre
Group
Missing Rows 7
Wilcoxon / Kruskal-Wallis Tests (Rank Sums)
Expected
Level Count Score Sum Score Score Mean (Mean-MeanO)/StdO
Post 92 4864.50 7084.00 52.875 -8.269
Pre 61 6916.50 4697.00 113.385 8.269
2-Sample Test, Normal Approximation
S Z Prob> IZI
6916.5 8.26872
1-way Test, ChiSquare Approximation
ChiSquare DF Prob>ChiSq
68.4025 1
Figure 21. Results of Wilcoxon Rank Sum test for TN at UCFRBA Station B3670000
The same analysis is conducted for TP (Figure 22) and turbidity (Figure 23). For TP, we cannot
reject the null hypothesis, and the p value is very high at 0.7209. In other words, there is no
statistically significant difference in TP following the plant upgrade. Mean TP concentrations
before and after the upgrade are 0.271 and 0.277 mg/L, respectively, while median TP
concentrations before and after the upgrade are 0.195 and 0.211 mg/L, respectively. For
turbidity, the p value is 0.0611, suggesting weak evidence of a change in turbidity (but not
significant at the recommended 95% confidence level); the decision to reject the null hypothesis
depends on what one selects as the critical p value. Mean turbidity readings before and after
the upgrade are 55.3 and 37.4 NTU, respectively, while median turbidities before and after the
upgrade are 30 and 19.4 NTU, respectively. The analyst should note the difference, but report
that the difference was not significant at the 95% confidence level.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 54

Oneway Analysis of TP By Group
2
1.5
•
� 1 =
•
0.5
0
Post Pre
Group
Missing Rows 9
Wilcoxon / Kruskal-Wallis Tests (Rank Sums)
Expected
Level Count Score Sum Score Score Mean (Mean-MeanO)/StdO
Post 83 6404.00 6308.00 77.1566 0.357
Pre 68 5072.00 5168.00 74.5882 -0.357
2-Sample Test, Normal Approximation
S Z Prob>IZI
5072-0.35725 0.7209
1-way Test, ChiSquare Approximation
ChiSquare DF Prob>ChiSq
0.1290 1 0.7195
Figure 22. Results of Wilcoxon Rank Sum test for TP at UCFRBA Station B3670000
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 55

Oneway Analysis of Turbidity By Group
2000
1000
'
700
500
400 •
300
200
' 100
70
50
40
~ 30
20
t
10
T
7
5
4
3
2
Post
Pre
Group
Missing Rows 7
Wilcoxon / Kruskal-Wallis Tests (Rank Sums)
Expected
_
Level Count Score Sum Score Score Mean (Mean -MeanO)/StdO
Post 92 6581.00 7084.00 71.5326
-1.873
Pre 61 5200.00 4697.00 85.2459
1.873
2-Sample Test, Normal Approximation
S Z Prob> IZI
5200 1.87252 0.0611
1-way Test, ChiSquare Approximation
ChiSquare DF Prob>ChiSq
3.5133 1 0.0609
Figure 23. Results of Wilcoxon Rank Sum test for Turbidity at UCFRBA Station B3670000
6.8 ARIMA Analysis
As noted at the conclusion of Section 4.5, this SOP recommends against using ARIMA for routine
trend analysis due to the complexity of the procedure, and because multiple good candidate
models can often be specified that predict different trends. However, ARIMA is also a valuable
tool for intensive special study analyses and an example ARIMA analysis is provided for
reference so that the steps in JMP can be documented.
In this example, NOx at Ellerbe Creek station EL1.9EC is considered. Data from ambient
monitoring only is used to prevent introduction of bias from special studies, which tend to
violate the ARIMA assumption of fixed -interval sampling, and may have different distributions
when alternate laboratories are used. First, a time series plot is provided (Figure 24).
Considerable variation is present, and there appear to be cycles present in the data that span
multiple years. One value below the detection limit of 0.1 mg/L is present, but is equal to the
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 56

minimum value of the entire series, so the value is set equal to the detection limit. Less than one
percent of the data are affected by detection limits.
Figure 24. Time Series of NOx at City of Durham Station EL1.9EC
A histogram of the data reveals some skew in the distribution (Figure 25). ARIMA is unlikely to
address unequal variance when orders of differencing are applied, so a log transformation of the
data should be considered. However, the skew of the log transformed data is higher than the
skew of the untransformed data, so the original data are used in the analysis.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 57

Figure 25. Histogram of of NOx at City of Durham Station EL1.9EC, 2004 — 2013
To start the ARIMA analysis in JMP, the Time Series Analysis window is opened by running
Analyze > Modeling > Time Series, specifying NOx as the input for Y, Time Series, and date as X,
Time ID. The Time Series windows opens, which shows the ACF and PACF, as described in
Section 4.5. The results suggest the presence of significant autocorrelation in the data.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 58

Time Series NOx
i
s
6
5
4
3
2
1
0
01/01/2004 01/01/2007 01/01/2010 01/01/2013
Date
Time Series Basic Diagnostics
Lag AutoCorr-.8-.6-.4-.2
0 ..2.4
.6 .8
0 1.0000
1 0.4127
2 0.3564
3 0.3103
4 0.2407
5 0.3460
6 0.2161
7 0.1672
8 0.1826
9 0.2086
10 0.1162
11 0.0649
12 0.0117
13 0.0315
14 0.0297
15-0.1225
16-0.1532
17-0.1854
18-0.1315
19-0.1767
20-0.2237
21-0.2142
22-0.1952
23-0.0973
24-0.0952
25-0.1218
Mean 2.4407692
Std 1.586448
N 117
Zero Mean ADF-3.268858
Single Mean ADF-6.830787
Trend ADF-6.855049
Ljung-Box Q
p-Value
Lag
0
Partial-.8-.6-.4-.2
1.0000 ;
0 .2
.4 .6 .8
20.4411
<.0001*
1
0.4127
35.8160
<.0001*
2
0.2242
47.5719
<.0001*
3
0.1315
54.7084
<.0001*
4
0.0404
69.5935
<.0001*
5
0.2118
75.4499
<.0001*
6
-0.0293 ;
78.9873
<.0001*
7
-0.0372
83.2441
<.0001*
8
0.0371
88.8528
<.0001*
9
0.1006
90.6086
<.0001*
10
-0.0989
91.1623
<.0001*
11
-0.0688 . • •
91.1803
<.0001*
12
-0.0626
91.3132
<.0001*
13
0.0085
91.4322
<.0001*
14
-0.0333
93.4794
<.0001*
15
-0.1756
96.7163
<.0001*
16
-0.1056
101.504
<.0001*
17
-0.0804
103.936
<.0001*
18
0.0088 ;
108.371
<.0001*
19
-0.0826
115.556
<.0001*
20
-0.0467
122.210
<.0001*
21
-0.0075
127.794
<.0001*
22
0.0099
129.197
<.0001*
23
0.1001
130.553
<.0001*
24
0.1051
132.798
<.0001*
25
0.0338
Figure 26. Time Series Analysis Window for NOx at City of Durham Station EL1.9EC
Various ARIMA models are run by selecting the down -pointing red triangle next to Time Series
NOx, and selecting one of the following options:
• ARIMA - for models of the form ARIMA(p,d,q). Leave Intercept checked to specify a
model with a constant.
• Seasonal ARIMA- for models of the form ARIMA(p,d,q)x(P,D,Q)12.
• ARIMA Model Group- multiple combinations of ARIMA and/or Seasonal ARIMA models
can be run by specifying starting and ending orders for each model term.
The Model Comparison window is shown below the initial ACF/PACF plots. Each ARIMA model
with Report checkbox selected is shown below the Model Comparison Window. Residuals can
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 59

be viewed for each model to check for stationarity. In this example, results for four models are
shown (Table 6).
Table 6. Statistics for Example ARIMA Models
Model
Variance
AIC
IMA(1,1)
2.05
411.8
ARMA(1,1)
1.95
413.5
ARIMA(1,0,1)x(0,1,1)12
2.06
400.9
ARIMA(0,1,1)x(0,1,1)12, no constant
2.27
391.6
The residuals of all of the models have ACFs with lags that immediately decay to low values,
although the IMA(1.1) model has lag coefficients of the lowest magnitude. The ARMA(1,1) has
the best variance among the models, while the ARIMA(0,1,1)x(0,1,1)12, with no constant has
the best AIC. Any of these models could be a successful candidate for consideration. However,
they all differ in terms of predicted trends.
Figure 27 shows the fit for each model along with the future forecast and confidence bounds on
the forecast. One visible difference is that the seasonal models include a seasonal pattern in
their forecast. Of greater concern is the difference in the direction and magnitude of
forecast/trend. The IMA(1,1) indicates an increase in NOx over time. Likewise,
ARIMA(1,0,1)x(0,1,1)12 predicts a subtle increase. The ARMA(1,1) model however predicts a
decrease, and the ARIMA(0,1,1)x(0,1,1)12, with no constant model indicates a decrease as well.
All the models have confidence bands that cover a range from increasing to decreasing trends.
Recall that the ARIMA forecasts include the effects of the persistence of random shocks as well
as any underlying deterministic trend. Indeed, only models that include the integrated
component (the third and fourth ARIMA models and not the IMA and ARMA models) include
any deterministic trend. In this case, the models without differencing (the IMA and ARMA
models) have lower variance and higher AIC; thus, there is little evidence for a deterministic
trend.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 60

IMA(1,1)
ARMA(1,1)
Forecast
Forecast
7
8-
A
a
6
•
6
•
q
• •
.+.yi �•�i,4
j
5
• •
• • •
2•
•�
• �� tb
c
2
' • •'•_•
0
•
• •
-2
i
01/01/2005 01101/2OD9 01/01/2013
01/01/2005 01/010009 01/01/2013
Date
Date
ARIMA(1,0,1)x(0,1,1)12
ARIMA(0,1,1)x(0,1,1)12 no constant
Forecast
Forecast
9
8-
fi
•
6 ••VV
a
•; •
•
•• • • P
•
2
yV
a
0
• • M '
-2
z
-4
01/01/2D05 01/01/2009 01/01/2013
01/01/2005 01/01f2009 01/01/2013
Date
Date
Figure 27. Forecast and Trend Predictions for Example ARIMA Models
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 61

(This page is intentionally left blank)
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 62

7.0 Quality Control
[City to insert its guidance on quality control]
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 63

(This page is intentionally left blank)
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 64

8.0 References
Box, G.E.P. and G.M. Jenkins. 1976. Time Series Analysis: Forecasting and Control, 2nd ed.
Holden -Day, San Francisco.
Box, G.E.P. and G.C. Tiao. 1975. Intervention analysis with applications to economic and
environmental problems. Journal of the American Statistical Association. 70:70-79.
Chandler, R.E. and E.M. Scott. 2011. Statistical Methods for Trend Detection and Analysis in the
Environmental Sciences. John Wiley and Sons, West Sussex, United Kingdom.
Cleveland, W. S. 1979.,Robust locally weighted regression and smoothing scatterplots.,Journal of
the American Statistical Association. 74:829-836.
Gilbert, R.O. 1987. Statistical Methods for Environmental Pollution Monitoring. Van Nostrand
Reinhold, New York.
Helsel, D. 2010. Much Ado About Next to Nothing: Incorporating Nondetects in Science. The
Annals of Occupational Hygiene. 54(3): 257-262.
Helsel, D.R., D.K. Mueller, and J.R. Slack. 2006. Computer Program for the Kendall Family of
Trend Tests. Scientific Investigations Report 2005-5275. U.S. Geological Survey. Reston, VA.
Hirsch, R.M., J.R. Slack. 1984. A nonparametric trend test for seasonal data with serial
dependence. Water Resources Research. 20:727-732.
Hirsch, R.M., J.R. Slack, and R.A Smith. 1982. Techniques of trend analysis for monthly water
quality data. Water Resources Research. 18:107-121.
Hirsch, R.M., R.B. Alexander, and R.A. Smith. 1991. Selection of methods for the detection and
estimation of trends in water quality. Water Resources Research. 27:803-813.
Hodges, J.L. and E.L. Lehmann. Estimates of location based on rank tests. Annals of
Mathematical Statistics. 34:598-611.
Kendall, M.G. 1975. Rank Correlation Methods. 4th ed. Charles Griffin, London.
Mann, H.B. 1945. Non -parametric test against trend. Econometrica. 13:245-259.
Mann, H.B, and D.R. Whitney. 1947. On a Test of Whether one of Two Random Variables is
Stochastically Larger than the Other. Annals of Mathematical Statistics 18(1): 50-60.
Nau, B. 2005. Forecasting, Decision 411.
http://people.duke.edu/—rnau/Decision4llCoursePage.htm. Accessed January 24, 2014.
Pindyck, R.S. and D.L. Rubinfeld. 1997. Econometric Models and Economic Forecasts, 41h ed.
Irwin/McGraw-Hill.
Sen, P.K. 1968. Estimates of the regression coefficient based on Kendall's tau. Journal of the
American Statistical Association. 63:1379-1389.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS MARCH 30, 2015
STANDARD OPERATING PROCEDURES PAGE 65

Tetra Tech. 2013. Lake B. Everett Jordan Watershed Model Report. Prepared for North Carolina
Nutrient Science Advisory Board, North Carolina Division of Water Resources, and Triangle J
Council of Governments.
CITY OF DURHAM WATER QUALITY TREND ANALYSIS
STANDARD OPERATING PROCEDURES
MARCH 30, 2015
PAGE 66