&EPA
United States
Environmental Protection
Agency
EPA 600/R-25/337 I April 2025 I www.epa.gov/research
Decadal Patterns of Climate Factors and Their
Influence on Vegetation Condition Across Five
Tillamook County (OR) Watersheds
Temporal trends in monthly Normalized Difference Vegetation Index
Tillamook basin, Oregon (1989 - 2019)
Univariate autoregression Multivariate autoregression
Office of Research and Development
Center For Public Health and Environmental Assessment
-------
EPA 600/R-25/337 I April 2025
Decadal Patterns of Climate Factors and Their
Influence on Vegetation Condition Across Five
Tillamook County (OR) Watersheds
by
Maliha S. Nash1, Matthew C. Harwell1, Megan VanFossen2, Daniel Rosenbaum2, and Theodore H.
DeWitt3
U.S. Environmental Protection Agency
Office of Research and Development
Center for Public Health and Environmental Assessment
Pacific Coastal Ecology Branch
Newport, Oregon 97365
1 U.S. Environmental Protection Agency, Newport, Oregon
2 U.S. Environmental Protection Agency, RTP, North Carolina
3 (retired) U.S. Environmental Protection Agency, Newport, Oregon
-------
EPA 600/R-25/337 I April 2025
Disclaimers/Notices
Disclaimers
This report has been subject to the U.S. Environmental Protection Agency, Office of Research and
Development peer and administrative review and approved for publication. Any mention of products,
services, agencies, organizations, or people does not imply an endorsement by the U.S. Government or
the U.S. Environmental Protection Agency.
Citation
The citation for this report is: Maliha S. Nash, Matthew C. Harwell, Megan VanFossen, Daniel
Rosenbaum, and Theodore H. DeWitt. 2025. Decadal Patterns of Climate Factors and Their Influence on
Vegetation Condition Across Five Tillamook County (OR) Watersheds. U.S. Environmental Protection
Agency, Office of Research and Development, Washington, D.C. EPA/600/R-25/337.
-------
EPA 600/R-25/337 I April 2025
Table of Contents
Disclaimers/Notices ii
Disclaimers ii
Citation ii
List of Tables iv
List of Figures iv
Acknowledgments viii
Abbreviations ix
Executive Summary 1
1 Introduction 2
2 Materials and Methods 5
2.1 Site Description 5
2.2 Data 7
2.3 Analytical Approach 8
2.3.1 Univariate Autoregression 8
2.3.2 Multivariate Autoregression 9
2.4 Analytical Interpretation 10
3 Analyses for Climatic Factors over 31 Years (1989-2019) in Tillamook Basin 12
3.1 Description 12
3.2 Seasonal Anomalies 18
3.3 Derived Climatic Factors: Apparent Temperature Index, Potential Evapotranspiration, and
Moisture Index 21
4 Trend Analyses for NDVI and related Climatic Factors over 31 years (1989-2019) in Tillamook Basin...27
4.1 Temporal Trend in NDVI and Climatic Factors over 31 years 27
4.2 Changes of NDVI including Climatic Factors 29
4.3 Climatic Factors Contribution 29
5 Cause of Greenness Change 31
5.1 Tree Mortality, Insect Infestation 31
5.2 Changes in Greenness in Restoration Areas 35
5.3 Changes in Greenness due to Climatic Factors 37
6 Major Findings and Conclusions 39
6.1 Major Findings 39
6.2 Conclusions 40
7 Quality Assurance 41
8 References 42
Appendix 1 46
Appendix 2 51
-------
EPA 600/R-25/337 I April 2025
List of Tables
Table 3-1. Monthly, overall, and seasonal summary statistics (minimum, mean, maximum, coefficient of
variation) for the climatic factors used in the study over 31 years for the Tillamook basin. Note: hPa =
hectopascal 12
Table 4-1. Percent of total areas with significant increasing or decreasing monthly trend in NDVI and
climatic factors from univariate autoregression model from 1989 to 2019 in the Tillamook basin, Oregon.
The last two rows are Percent of total areas with significant increasing or decreasing trend in NDVI and
from univariate and multivariate autoregression models from 1989 to 2019 in the Tillamook basin,
Oregon 27
Table 4-2. Frequency of outcomes in comparison of NDVI trends between univariate and multivariate
analyses for 3,291,543 30 m x 30 m pixels in the Tillamook basin, Oregon. Outcomes A and C (bold)
denote changes due to local factors; Outcome B denotes changes due to climatic factors, and Outcome D
denotes no change in greenness 30
Table 5-1. Univariate trend for climatic factors for Sample Locations (SL) within a tree damaged sample
(pixel: 30 m x 30 m). Tree damage was observed in sample (year): Sample 1 (1996); Sample 2 (1991,
1993); Sample 3 (1991, 1992, 2006, 2016); Sample 4 (2019); Sample 5 (1991, 1992, 2004, 2006, 2015);
Sample 6 (1993); Sample 7 (1991, 2006, 2013, 2015); Sample 8 (2013). Vertical shaded columns denote
that tree damaged was observed in more than one year 35
Table App-1. Summarized totals for geographically corresponding unfiltered and filtered Landsat scenes
for each Landsat platform (with GEE Catalog identifier) 52
List of Figures
Figure 2-1. Landcover in 2001 and 2019 for the Tillamook basin (A, B), watersheds and streams in the
basin (C), and (D) elevation for 10 m pixels, shaded areas denote elevation from O-m(yellow) to 1125-m
(dark blue). Landcover classes were lumped showing the most common areas. White shaded areas in (A)
and (B) denotes the clear-cut areas. Light cyan shaded areas in (C) denote Tillamook Bay estuaries
(woody wetland and emergent herbaceous wetlands) 6
Figure 2-2. Graphical abstract for our approach in comparing NDVI trend between univariate and
multivariate regressions. The four outcomes are: Outcome A (local change): slope of time is significant
from both univariate and multivariate models. Outcome B (climate signal): slope of time is significant in
univariate model but not in multivariate model. Outcome C (local change): slope of time is not significant
in univariate model, but it is significant in multivariate model. Outcome D (no change): slope of time is
not significant in both univariate and multivariate models 11
Figure 3-1. Monthly precipitation within the Tillamook basin over 30 years (1989-2019; n=372 months).
(A) average precipitation (mm/month), (B) coefficient of variation [CV] of monthly precipitation. Average
wet season total precipitation (C) (mm/season) reveals that wet average precipitation values are higher
in the middle of the basin and lower within the western and eastern bounds. (D) coefficient of variation
[CV] of wet season precipitation. Average dry season total precipitation (E) (mm/season). (E) (coefficient
of variation [CV] of seasonal precipitation) shows that variability <50% for dry precipitation; CV is higher
at the eastern side, then decreases gradually toward the western side of the basin 14
iv
-------
EPA 600/R-25/337 I April 2025
Figure 3-2. The temporal trend of monthly precipitation in 30 m x 30 m pixels over 30 years (1989 -
2019). Map A (trend values mm/month of monthly precipitation) shows that trend values decreased
within the basin and ranged from -0.052 to -0.251 mm/month. Map B (direction and probability of
monthly precipitation) shows that monthly precipitation decreased significantly (p < 0.05) in the middle
of the basin. Rainfall station, South Fork (GHCND:USROOOOOSFK), is in Map B within the area of non-
significant precipitation trend. Slope direction and slope probability of precipitation from the rainfall
station is same as of that from Map B (non-significant grey are in the map trend from PRISM data).
Seasonal precipitation trend values within the Tillamook basin over 30 years (n=30 years). Map C (slope
value of wet precipitation, mm/season) shows that the magnitude of wet trends on the eastern and
western part of the watershed are less than those in the middle of the watershed in both seasons. Map D
(slope value of dry season precipitation, mm/season) reveals that dry season precipitation decreased
with years but with lower rate than that of wet season, all are within the low class of wet season (< 5
mm/season). Prediction in legend biplot is predicted NDVI 15
Figure 3-3. Monthly precipitation (PPT) and potential evapotranspiration (PET) for 2006, 2007, and 2008.
Trees were removed after 2006. The location of this pixel is marked by the cyan symbol in Figure 6. Pixel
size is 30- x 30-m 16
Figure 3-4. 2008 National Land Cover Database (NLCD; h ftps://www.mrlc.gov/data) of Tillamook basin;
the cyan colored point is within a clearcut area. Landscape type in areas where the marked point (Point
25) is located was evergreen in 2001 -2006, grass (71) in 2008-2016, then shrub and scrub (52) in 2019.
Biplot of NDVI i/s time (inset) shows a dip around 2007. The overall decrease in NDVI over 31 years was
significant 17
Figure 3-5. Monthly precipitation (PPT) and potential evapotranspiration (PET) for 2006, 2007, and 2008.
Trees were intact in 2006, 2007, and 2009. The location of this pixel is marked on the Google Earth map
in Figure App-4. Pixel size is 30- x 30-m 18
Figure 3-6. Seasonal precipitation (mm/season) over 30-years for the selected pixel. The horizontal line at
0 in both figures represents the seasonal precipitation 31-year average. (A) wet season anomaly
precipitation and (B) dry seasonal precipitation anomaly. Seasonal precipitation anomaly > 0 values for
wet and dry season. (C) Wet season precipitation anomaly values and its trend (y = -7.2745x + 15167, R2
= 0.0194). (D) Dry season precipitation anomaly values and its trend (y = 5.7265x -11335, R2 = 0.19).... 19
Figure 3-7. Number of years with precipitation anomalies. (A) Wet season number of years > 30-year
average (2,653.9 mm/season). (B) Dry season number of years > 30-year average (391.1 mm/season). (C)
Wet season number of years < 30-year average (2,653.9 mm/season), and (D) Dry season number of
years < 30-year average (391.1 mm/season) 20
Figure 3-8. Distribution of wet and dry seasonal anomaly trends. (A) Wet anomaly > 30-year wet season
average. (B) Dry anomaly > 30-year dry season average. The darker shaded areas denote higher trend
values of seasonal precipitation (amount) anomalies. Estuaries experiences more precipitation within the
dry season but experience no change from the wet season. Blank areas indicate no trend
evaluation 21
Figure 3-9. PET annual average (A), precipitation annual average (B), and annual average of moisture
index (C). The darker shaded areas indicate higher average values. Higher annual PET are coinciding with
v
-------
EPA 600/R-25/337 I April 2025
that of less precipitation. Moisture index values ranged from 0.37 to 0.76 which falls in the wet to very
wet scales per Feddema (2005) wetness classification. The distribution of moisture index in the basin
corresponds with that of precipitation pattern. Pixel size is 30 m x 30 m 22
Figure 3-10. Spatial distribution of average (n=31 years) precipitation (PPT, mm/month), potential
evapotranspiration (PET, mm/month), and moisture index (Ml, unitless) per month. Scale of Ml
(Feddema, 2005; 1.00 =>MI >= 0.66 very wet, 0.33 <= Ml < 0.66 very wet, 0<= Ml < 0.33 moist, -0.33 <=
Ml <0 dry, -0.66 <= Ml < -0.33 semi-arid, -1.00 <=MI < -0.66 arid). Pixel size is 800x800 m 26
Figure 4-1. Univariate autoregression model results: (A) trend values range between -0.0024 to + 0.0022,
(B) significant temporal trend for the monthly NDVI from 1989 through 2019 in the Tillamook basin,
Oregon. Green pixels indicate significant (p < 0.05) increase in greenness (i.e., NDVI); red pixels indicate
significant decrease in greenness 27
Figure 4-2. From univariate autoregression model: pixels (30 m x 30 m) with significant temporal trend
(1989-2019) for monthly: (A) precipitation (red indicates significant [p < 0.05] decrease [drier]); (B)
maximum temperature (blue indicates significant decrease [cooler], red indicates significant increase
[warmer]); (C) minimum temperature (red indicates significant increase [warmer night temperature]);
and (D) maximum vapor pressure deficit (red indicates significant increase and blue indicates significant
decrease). In all maps, grey indicates no significant change 28
Figure 4-3. Significant temporal trend for the monthly NDVI from 1989 through 2019 in the Tillamook
basin, Oregon. (A) univariate autoregression model and (B) multivariate autoregression model results:
pixels (30 m x 30 m) with green pixels indicate significant (p < 0.05) increase in greenness (i.e., NDVI); red
pixels indicate significant decrease in greenness 29
Figure 4-4. The spatial distribution (left) of the four outcomes (right) from comparison of the NDVI trend
significance between univariate and multivariate autoregression for each pixel. Definition of the four
outcomes is in the contingence table (right). Both Outcomes A & C reflect changes due local factors,
Outcome B reflects a climate factor, and Outcome D reflects no change 30
Figure 5-1. New incidence of tree mortality for (A) 1991, (B) 2006, (C) 2013, (D) 1990-2019 cumulative
tree mortality. Map with color polygons over the years in Figure App-4 32
Figure 5-2. A significant decreasing trend (slope = -0.00108 p < O.OOOl) between maximum NDVI and tree
mortality (DAM1) from 1989 through 2019 in the Tillamook basin 33
Figure 5-3. Locations of tree mortality samples. Marked polygons include tree damage from one year to
several. Tree mortality observed in Sample 1:1996; Sample 2:1991, 1993; Sample 3:1991, 1992, 2006,
2016; Sample 4:2019; Sample 5:1991, 1992, 2004, 2006, 2015; Sample 6:1993; Sample 7:1991, 2006,
2013, 2015; Sample 8:2013. The inset biplot shows the NDVI temporal behavior in a pixel (row=1692
col=210) within Sample 8, with abrupt change around 2008-2009 34
Figure 5-4. Top-right, location of the restoration polygon within the Tillamook Basin; top-left, green,
orange, and grey polygons denote native dominant vegetation, non-dominant vegetation, and bare
ground, respectively, as described by Janousek etal. (2021). Bottom, selected restoration polygons from
Janousek et al. (2021) overlay two of the four outcomes (A & D). Outcome D sample is within the
vi
-------
EPA 600/R-25/337 I April 2025
northern zone of the SFC and Outcome A is in the southern zone of the SFC. Elevation gradient sorted as
Outcome sample D is at higher elevation than sample A. The locations of the outcome samples on Google
Earth and the biplot ofNDVI i/s months are in Figures App-3 and App-4for Outcome D sample. Biplot of
NDVI i/s months for Outcome sample A is in Figure App-5 36
Figure 5-5. Outcomes related to climatic factors (top, green denotes local changes, dark red denotes
changes related to climatic factors) 37
Figure 5-6. Historical NLCDfor areas around sample raster (Cyan closed circle). This sample remains
within a tree intact area in 1994, 2003, 2005, 2011, 2016, 2019, and 2024, whereas tree removal to the
west of the sample location took place between 2005 and 2011 38
Figure App-1. Landcover/landuse from 2001 to 2019 in Tillamook basin, Oregon 46
Figure App-2. Monthly significant trend within the Tillamook basin over 30 years (1989-2019; n=372
months) for average temperature (°C) (A), dewpoint temperature (°C) (B) and minimum vapor pressure
deficit (hPa) (C). Red shaded areas denote significant increase, shaded green areas denote significant
decrease 47
Figure App-3. Google Earth map showing location of sample for an Outcome D in 2005 (A) and 2011 (B).
Trend of monthly NDVI (C) for the sample in (A) and (B) decreased non significantly (p=0.057) and
remained not significantly decreasing when climatic factors were included (P=0.07) 50
Figure App-4. Colored polygons show the variability in tree mortality outbreak timing from 1989 to 2019.
Tree mortality ranged from 1 (1989 and 1990) to 100 (2008) number of dead trees per acre. Other years
with high tree mortality, number of dead trees per acres for 1992 and 2018 were 62 and 41, respectively.
49
Figure App-5. Google Earth map for Outcome A sample location and its NDVI i/s months 50
VII
-------
EPA 600/R-25/337 I April 2025
Acknowledgments
We are very grateful to several people who contributed crucial information, ideas, support,
encouragement, or constructive review that helped improve the quality of this report. These people
include Dr. Patti Meeks, Dr. James Wickham, and Dr. Andrew Pilant. Valuable external technical peer
reviews were contributed by Drs. Brady Couvillion, Heather Tollerud, and Melanie Vanderhoof.
viii
-------
EPA 600/R-25/337 I April 2025
Abbreviations
ACF
Autocorrelation Function
AT
Apparent Temperature Index
AVHRR
Advanced Very High-Resolution Radiometer
BQA
Band Quality Assessment
°Cday
Degree days in Celsius
CV
Coefficient of Variation
DW
Durbin-Watson statistic
EE
Earth Engine
EPSG
European Petroleum Survey Group
GEE
Google Earth Engine
hPA
Hectopascal
Ml
Moisture Index
NDVI
Normalized Difference Vegetation Index
NLCD
National Land Cover Database
NEP
National Estuary Program
OLI
Operational Land Imager
PET
Potential Evapotranspiration
PPT
Monthly Precipitation
PRISM
Parameter-elevation Regressions on Independent Slopes Model
SR
Surface Reflectance
TM
Thematic Mapper
Tmax
Monthly Maximum Temperature
Tmin
Monthly Minimum Temperature
TOA
Top-of-Atmosphere
VPD
Vapor Pressure Deficit
ix
-------
EPA 600/R-25/337 I April 2025
Executive Summary
The purpose of this report is a long-term monitoring of an environmental index, which is an essential
element in natural resources management. Conservation of natural resources, for example, requires a
well-established plan that includes a monitoring program to determine outcomes, and verify and
validate these outcomes. These analyses can provide researchers and environmental decision-makers
with early warning signals for widespread general trends as well as a means to identify specific areas
where land conditions are degrading or improving. In this report we describe long-term changes in
greenness and climatic factors for the Tillamook basin from 1989 to 2019. Changes in greenness across
terrestrial landscapes can result from human activities or natural processes. The rate of losing or gaining
greenness may occur over short or long periods of time depending on the individual activity or
combined activities. Long-term data sets of greenness typically cover large landscapes, and when paired
with weather and land-use data, provide an opportunity to quantify spatial and temporal changes in
vegetation, and their rates. This analysis of temporal trends of greenness is conducted to quantify
changes, and the direction of these changes, in the Tillamook basin using statistical models. Rapid
changes in greenness were found to indicate the occurrence of local factors such as fire, agriculture, and
removal of greenspace, whereas gradual changes in greenness can be related to changes in climatic
factors directly, or indirectly such as tree mortality from disease.
The report describes the temporal trend in selected climatic factors and greenness and their correlation
using time series analyses. Results summarize changes in greenness as related to local or climatic factor
changes. Within the Tillamook basin, several areas were examined via Google Earth and leveraging local
expert knowledge to validate the cause of change in greenness. The report describes and characterizes
climatic factors (minimum, maximum, mean, and dewpoint temperatures, and minimum and maximum
vapor pressure deficit, and derived climatic factors (potential evapotranspiration, moisture index, and
apparent temperature)) over 31 years. The report also describes seasonal precipitation (wet and dry
seasons) and seasonal anomalies (seasonal total precipitation > seasonal average) that were analyzed to
characterize trends and other statistics. The statistical analyses were applied in each pixel and results
were mapped to visualize changes in the Tillamook basin.
1
-------
EPA 600/R-25/337 I April 2025
1 Introduction
Abrupt changes in long-term trends of climatic factors can represent climatic signals that have
consequences directly or indirectly on both the environment and human health. Heat, drought, wildfire,
and floods are examples of climatic signals. Nighttime temperatures are getting warmer and warmer
temperatures are responsible for increases in drought condition (Weiss et al., 2009; Nash and
Christensen, 2017; Williams et al., 2020; Bonfils et al., 2020) and declining snow cover for the last 90
years (Najafi et al., 2016). Higher temperatures increase evapotranspiration and decrease available soil
water content causing wilting and drying of vegetation that can promote wildfire (Abatzoglou and
Williams, 2016).
Climatic signals can be detected by analyzing long-term time series data. Time series data can be from
observational data (e.g., rain gauges), model results, remote sensing data, and other sources. Coupling
climatic factors with that of human and environmental data in advanced statistical models can
determine relationships and changes. Monitoring vegetation change (e.g., greenness) over time at the
same location enhances understanding of status, condition, and trends by identifying changes in the
shape, size, and condition of vegetation areas. Remote sensing satellites provide a greenness index
(normalized difference vegetation index; NDVI) over space and time, enabling researchers to follow
vegetation changes in diverse areas. Accurate identification of plant type and phonologies of grassland
and shrubland were distinguished from the time series analyses of greenness index in the Jornada
Desert, New Mexico (Peters et al., 1997). The NDVI has been used to quantify variation in seasonal net
carbon fixation and net primary productivity and its relationship to anthropogenic sources on different
scales (Chhabra and Dadhwal, 2004). Moreover, NDVI has been used to quantify monetary losses in crop
yields due to abrupt stressors (Deininger et al., 2022). The NDVI can be used over diverse areas to
quantify the available forage and estimate grazing days left for each field (Flynn, 2006) for crop
insurance (Turvey and Mclaurin, 2012). Lybbert et al. (2011) evaluated argan tree NDVI changes in
Morocco using commune-level data on educational enrollment over the period from 1981 to 2009. The
LADAMER program (2003) mapped results of temporal trends of NDVI for decision-makers and
policymakers in Europe by utilizing trend maps to locate changes in vegetation cover (positive and
negative) on different scales. Special attention and further investigation were given to areas with
negative trends to determine if the cause related to anthropogenic activities. From a monitoring study in
an urban environment, Miller et al. (2022) reported the type of correlation between NDVI, temperature,
and drought and their effect on different vegetation types (trees and grass covers) and those
correlations were different between season and vegetation types.
Long-term, spatiotemporal analyses of climatic factors can have a significant role in studies related to
environmental changes, human health, availability of ecosystem services, and others. Many studies
demonstrated the significant association of selected climatic factors with changes in the environment
(Bonfils et al. 2020), shifts in health outcomes (Perry et al., 2011), and fog formation (Hiatt et al., 2012).
The long-term data sets used in this study can be used to examine trends in temperature and
precipitation to answer environmental management questions such as, "What is the influence of
drought, heat waves, and potential evapotranspiration on potential shifts in water storage, hydrological
flows, biogeochemical processes, or biological shifts on water-dependent species in the Tillamook
basin?" Vanderhoof et al. (2018) utilized a summarization of the 25-year climate and derivative datasets
from Nash and Christensen (2017) in a drought resilience study to explain changes in water storage as
measured by inundation maps.
2
-------
EPA 600/R-25/337 I April 2025
Changes in greenness across terrestrial landscapes can result from each or combined effects of human
activities and natural processes. The rate of losing or gaining greenness may occur over short or long
periods of time depending on the type of individual or combined activities. Long-term data sets of
greenness typically cover large landscapes, and when paired with weather and land-use data, provide an
opportunity to quantify spatial and temporal changes in vegetation, and their rates. The analyses of
temporal trends of NDVI data using statistical models quantifies changes and their direction in a given
watershed or basin. Rapid changes in NDVI indicate the presence of local factors such fire, agriculture, or
removal of greenspace, whereas gradual changes in NDVI can be related to changes in climatic factors
directly, or indirectly, such as tree mortality from disease. Drought and biological resistance due to
changes in climatic factors reduce trees resistance to insect infestation (Bernal et al., 2023).
Teasing apart climate-related changes (i.e., gradual shifts in NDVI) from other stressors (e.g., direct and
rapid changes due to fire, land use change, etc.) is challenging as changes in NDVI can be confounded by
time-dependent patterns (e.g., seasonal leafiness, vegetation succession) and variation associated with
climatic factors. The use of time series regression models that account for serial autocorrelation
improves our ability to identify significant changes of greenness with time. In some studies, the Durbin-
Watson (DW) statistic was evaluated to test for the presence of serial correlation (e.g., Ji and Peters,
2003), the non-parametric Mann-Kendall test was used to test for the existence of a monotonic trend,
and Theil-Sen statistic was used to account for the effects of serial correlation in detection of trends, but
has not been used as widely as DW and Mann-Kendall (Alashan, 2019). In contrast to the above
methods, autoregression regression model can be univariate or multiple regression where the trend of
response variable is adjusted for presence of covariance(s) in the model and the significant value of the
slope is adjusted for serial correlation that a rise from temporal variability.
The Tillamook basin is within the Oregon coastal and estuaries where the ecosystem is experiencing
changes in greenness, water quality, biota populations, and many other system qualities that are
thought to be attributable (directly or indirectly) to climate change (Ruggiero et al., 2010; TBNEP 2020).
An average from 13 Global Climate Models (GCMs) were scaled down to local scales (fine scale) used to
determine trends in climatic factors over time. Model projections average temperature for summers will
increase by 5-8°F and winter 4-7°F by 2085 in Tillamook estuaries. Past and future projection of climatic
factors were presented in Koopman (2018). Algal blooms, lower oxygen levels, and others are associated
with higher temperature which has direct effect on natural resources. The objective of this study is to
apply time series analyses autoregression models to analyze long-term changes in NDVI across the
Tillamook basin at a 30-m resolution, incorporating other data sources to examine multiple climatic
factors. In contrast to using an ordinary least square regression model to identify significant changes in
NDVI, we used an alternative method as autoregression models to improve our ability to identify trends
in NDVI with time. The trend in NDVI combining novel univariate (NDVI = time) and multivariate
(NDVI=time + climatic factors) autoregression models allow us to relate greenness trend to specific
climatic factors or non-climatic factors across the basin. The analyses were done per pixel following
previous approaches where temporal NDVI data for the state of New Mexico, USA (Nash et al., 2014),
and the continental United States (Nash et al., 2017) were used. The purpose of this report is to analyze
changes in greenness in the Tillamook basin. We use a 31-year time series of NDVI to analyze long-term
changes in greenness using single (time only) and multivariate autoregressive models to account for
spatial variation in the relative importance in abrupt (land conversion, fire) and longer-term changes in
greenness.
3
-------
EPA 600/R-25/337 I April 2025
Our research questions can be summarized as:
1. Are temperatures, vapor pressure, and derived factors (potential evapotranspiration, moisture
index and apparent temperature) over the 31 years increased or decreased significantly from
1989 to 2019 in the Tillamook Basin?
2. Is NDVI decreased or increased significantly from 1989 to 2019 in the Tillamook Basin?
3. Is seasonal precipitation (wet and dry seasons) and seasonal anomalies (seasonal total
precipitation > 31 seasonal average) trend over the 31 years increased or decreased significantly
from 1989 to 2019?
4. How are NDVI trends different from univariate and multivariate models? and
5. Specify locations of pixels with changes related to climatic factors or not.
Our presentation of the long-term datasets and trends in precipitation, temperatures, and vapor
pressure deficit can be used to study the influence of these climatic variable trends on greenness, water
storage, biological shifts, water-dependent species, and others. Our intent was to make the
methodology and analyses available for research efforts such as these, to better understand the trends
of NDVI and climatic factors and their interactions in hydrological, societal, biological, and ecological
systems within the Tillamook Basin.
4
-------
EPA 600/R-25/337 I April 2025
2 Materials and Methods
As is the case elsewhere, ecosystems in the Tillamook basin are experiencing the stress of climatic
changes on water quality, trees, and marine species that have an economic value for people (Ruggiero
et al., 2010). This section describes the methodology, data, and analyses available to support other
research and environmental management efforts needed to understand the trends of precipitation,
temperature, and their interactions in hydrological, societal, biological, and ecological systems in the
Tillamook basin. Results from this study will help in assessment, predicting, and defining management
actions in conserving the natural ecosystems and socioeconomic well-being of all residents. Further, the
approaches used here are fully transferrable to other systems (e.g., the other 27 National Estuary
Programs; NEP), spatial scales, and decision contexts.
2.1 Site Description
The Tillamook basin (south Tillamook County, Oregon, USA) consists of five watersheds covering nearly
1500 km2 (Miami, Kilchis, Wilson, Trask, and Tillamook) with the main outlet of the basin draining to the
Pacific Ocean at Tillamook Bay (Figure 2-1C). The elevation varies from 0 feet (0 m) within the estuaries
to more than 2,000 feet (610 m) in the eastern part of the basin (Figure 2-1D). From NLCD GIS maps,
evergreen and mixed forests are the dominant landcover (~70%) in the basin with high density in the
eastern part that decreases gradually toward the west (Figures 2-1A and 2-1B). Fifteen percent of the
basin is developed, shrub/scrub, grassland/herbaceous and pasture. Landcover and land use data in
Figures 2-1A and 2-1B shows the changes for time period of between 2001 and 2019. Evergreen and
mixed forest decreased by 4% (from 73.55% to 69.61 %) from 2001 to 2019. Shrubland, grassland and
developed land (excluding developed mixed) increased by 2.26%, 1.73%, and 0.11, respectively. Woody
wetland, emergent herbaceous wetland did not change from 2001 to 2019. Pasture/hay decreased by
0.05% and water increased by 0.05%. For more temporal sequences to visualize changes in
landcover/land-use, Figure App-1 presents the landcover/land-use from 2001 through 2019, every 2 or 3
years.
5
-------
EPA 600/R-25/337 I April 2025
A:2001
B:2019
11, water
¦ 21, 22, 23, 24, Devloped
¦I 42, 43, Mixed & Evergreen
81, Pasture/Hay
82, Cultivated Crops
90, Woody Wetlands
95, Emergent Herbaceous Wetlands
C: Watershed in Tillamook basin
D: Elevation
I Kilometers
0 5 10 20 30 40
Figure 2-1. Landcover in 2001 and 2019 for the Tillamook basin (A, B), watersheds and streams in the basin (C), and (D)
elevation for 10 m pixels, shaded areas denote elevation from O-m(yellow) to 1125-m (dark blue). Landcover classes were
lumped showing the most common areas. White shaded areas in (A) and (B) denotes the clear-cut areas. Light cyan shaded
areas in (C) denote Tillamook Bay estuaries (woody wetland and emergent herbaceous wetlands).
6
-------
EPA 600/R-25/337 I April 2025
2.2 Data
A suite of remote sensing and environmental data were analyzed, including: NDVI; monthly averages of
minimum, maximum, average and dewpoint temperatures; precipitation; minimum and maximum vapor
pressure deficit; and tree mortality. The monthly NDVI data were obtained from Google Earth Engine
(GEE) (period of record: January 1989-December 2019) which was derived from LANDSAT data at a
resolution of 30 m x 30 m (a total of 3,291,543 pixels). The 31 years is sufficient to uncover periods of
NDVI gradual or abrupt change in the Tillamook Basin, this period of records agrees with that of
Anyamba and Tucker (2005), and Tucker and Anyamba (2011). Detailed description of the NDVI GEE
data is presented in Appendix 2.
This section concentrates on the potential of past, present, and future uses for long-term climate factor
trends, and trends of derived climatic factors and indices in the Tillamook basin. Analyses of all factors
over the 31 years (1989-2019) show the magnitude and direction of significant changes on different
time scales (monthly, annual average, and seasonal). Additionally, annual averages, coefficient of
variation (CV), minimum, and maximum values are also presented to show the differences in patterns of
variability and extreme values. The climate factors considered are: minimum temperature; maximum
temperature; dew point temperature; minimum and maximum vapor pressure deficit; and precipitation.
Additional, derived climatic factors were considered, including potential evapotranspiration (PET) and
two climatic indices (a moisture index (Ml) and an apparent temperatures (AT) index, often referred to
as a heat index) (Smoyer-Tomic and Rainham, 2001). Potential Evapotranspiration (PET) represents the
amount of water loss via plant and soil transpiration. PET is estimated from temperature, water
saturated vapor pressure and daylight length which varies with season and spatial locations. Ml is a
composite of precipitation, day length, average temperature and saturated vapor pressure and
represents the relative wetness (+ Ml values) and dryness (- Ml values) of an area. Apparent heat (AT) is
a heat index and has a relation to human health, derived by Smoyer-Tomic and Rainham (2001) (also
known as the index of heat; Perry et al., 2011) combines daily temperature and humidity. Derivation of
PET, AT, and Ml are described in detail in Nash and Christensen (2017).
The climatic parameters for the Tillamook basin analyses were derived from the PRISM (Parameter-
elevation Regressions on Independent Slopes Model) dataset (PRISM Climate Group;
http://prism.oregonstate.edu) at a scale of 800 m x 800 m. The 800 m x 800 m PRISM pixels were
downscaled into 30 m x 30 m using inverse distance weighting to match the spatial resolution of the 30
m Landsat-derived NDVI Google Earth Engine data. We chose to use a finer-resolution, 30-m grid scale
as it matched the resolution used by the VELMA (Visualizing Ecosystem Land Management Assessments)
eco-hydrological modeling platform that has been used extensively throughout the Pacific Northwest
(https://www.epa.gov/water-research/visualizing-ecosvstem-land-management-assessments-velma-
model) and is at a resolution meaningful for one of our target report audiences (the Tillamook Estuaries
Partnership). In previous study (Nash et al., 2014) we downscaled 4-km2 PRISM climatic variable to
match 1-km2 NDVI. Downscaling was also applied by Thorne et al. (2012).
The climatic factors provided by the PRISM data included monthly averages of daily minimum,
maximum, average and dewpoint temperatures (°C), minimum and maximum vapor pressure deficit
(hPa), and monthly total precipitation (mm). Previous research has shown that greenness is strongly
associated with temperature and precipitation (Wang et al., 2003; Twumasi et al. 2011). Vapor pressure
7
-------
EPA 600/R-25/337 I April 2025
deficit (VPD) is a measure of atmospheric dryness that influences vegetation (Konings et al., 2017),
where higher VPD results in higher moisture loss from plants. Konings et al. (2017) also indicated that
plants are more influenced by atmospheric dryness (i.e., VPD) than precipitation. Higher temperature
increases VPD. All analyses were performed for each of the 30 m x 30 m pixels for the Tillamook basin.
Maps of spatial distribution of temporal changes of climatic factors and derived climatic factors (e.g.,
PET) are presented below. Additional analyses were performed on derived climatic factors; PET, Ml, and
AT for each 800 m x 800 m pixel due to unclear patterns of spatial distribution of Ml and AT.
Tree mortality data were from the US Forest Service survey (USDA Forest Service, 2016) where yearly
polygons of infested areas (data collected early July through September) were mapped from 1989 to
2019 in the Tillamook basin. The mortality variable 'DAM1' was used to represent the number of dead
trees per acre. Tree mortality polygons were rasterized to 30 m x 30 m pixels and synchronized with that
of NDVI. Pixels with yearly maximum NDVI and tree mortality data were extracted to examine the
association between NDVI and tree mortality.
2.3 Analytical Approach
Two statistical time series model were conducted for each of the 3,291,543 pixels in the Tillamook basin.
First, univariate autoregression of NDVI versus time was conducted to quantify the temporal trend
(slope) for NDVI, and for each of the climatic parameters (e.g., precipitation versus time). The trend
direction for each significant pixel was then mapped to identify geographic patterns of significance and
trend direction. Second, a multivariate autoregression of NDVI versus time and climatic parameters was
conducted to reveal the combined effects and relative contributions of climatic factors to significant
NDVI trends. Time-series regression (autoregression) was used in both analyses because errors in
temporal data may be serially dependent, and if dependency exists, the standard error of the estimate
(e.g., slope) would be inflated. The autoregressive error model included a backstep function with up to
specified number of lags to account for serial dependence in errors. Autocorrelation function (ACF) and
partial autocorrelation function for residuals were checked for no significant correlation. A detailed
description of the model is described below. SAS/Stat software was used for all analyses.
Overall summary statistics were first calculated for each 30 m x 30 m pixel in the basin, followed by a
trend analysis. Additionally, the trend analyses for average temperature, dewpoint temperature, and
minimum vapor pressure that were not included in the NDVI climate factors section, are presented in
Figure App-2. For precipitation, additional analyses were conducted on a monthly and seasonal basis,
and for wet (October through April) and dry (May through September) seasons. Trend analyses applying
univariate autoregression were conducted for each of the climatic parameters (e.g., precipitation versus
time). The trend direction for each significant pixel was then spatially mapped to identify any geographic
patterns of significance and trend direction. Detailed model description is provided in the Appendix.
2.3.1 Univariate Autoregression
This analysis addressed trends in monthly NDVI and the individual climatic factors over the 30-year
period. For each 30 m x 30 m pixel, the autoregression model (Proc Autoreg; SAS/ETS, 2016) with
stepwise selection for the significant autoregressive error was fitted to the observed values to define the
direction and p-value for the slope as (Equation 2-1):
8
-------
EPA 600/R-25/337 I April 2025
Yt = 0o + 61* Time + pit (2-1)
and
Ht = T.t=iPil*t-i + £t (2-2)
Yt is and individual time series variable: monthly NDVI and any of the monthly climatic factors (Time= 1 -
372 months). The fitted autoregression models for the observed variable Yt, is the structural part plus
the autoregressive error (m). The structural part, which is the same as that of an ordinary least square
regression model (OLS; 0O + 0i *Time). The autoregressive error model (Equation 2-2) will account for
serial autocorrelation. The term Ht=i Pi flt-i 'n Equation 2-2 is the summation of the significant
autoregressive parameter (p) times lagged error(s), and k is the order of significant lags in the model.
The slope (0i) quantifies the rate and direction of change for each variable over 31 years. A significance
level of p < 0.05 was used to test whether the slope is differed from zero.
2.3.2 Multivariate Autoregression
The multivariate model for each pixel was:
NDVIt = /?0 + [>2Tmint + [>3Tmaxt + f3AVPDmaxt + [>5Timet+ fit (2-3)
and
Ht = T.t=iPil*t-i + £t (2-4)
where Ptis precipitation at month t, Tmin is minimum temperature, Tmax is maximum temperature, VDPmax
is maximum vapor pressure deficit, and e is the error term. The right side of Equation 2-3 includes the
autoregression error (nt) and the structure term (the remainder of the model), and by summing both
these terms yields the predicted value. The estimates (/?/s) for each factor quantify the magnitude and
direction of the relationship between NDVI and each factor over the 30-year period. The coefficient of
Time (f>5) is the temporal trend of NDVI after accounting for climatic factors. A comparison between
time coefficients of univariate (0i) with that from multivariate (|3s) models will discriminant between
whether significant NDVI trend is associated with climatic factors (and hence climate change) or local
factors (e.g., insect infestation). Four possible outcomes and inferences from the comparison are
summarized in two-way contingence table (Figure 2-2). Blocks of pixels that had either significant
increase or decrease (p < 0.05) in greenness were chosen to aid interpretation of results by drawing on
available literature, consultation with local experts, and Google Earth™. Average temperature, average
dewpoint temperatures, and minimum vapor deficit were not used in Equation 2-3, but their spatial
trend distribution may help explaining changes in NDVI.
9
-------
EPA 600/R-25/337 I April 2025
2.4 Analytical Interpretation
Changes in NDVI trend that can be linked directly to factors such as fire or landcover change, or
indirectly linked to factors such as climatic factors, can be assessed by comparing the results of the
univariate and multivariate analyses for each pixel. The comparison between the univariate and
multivariate models was summarized in a 2 x 2 contingency table (Figure 2-2). A significant NDVI trend is
indicated by a significant time coefficient, i.e., 0i in the univariate analysis (Appendix 1; Equation 2-1)
and (B5 in the multivariate analysis (Appendix 1; Equation 2-3). The four possible outcomes related to
specific changes are explained in Figure 2. Outcome A is a trend that can be interpreted as driven by
local changes, and occurs where the NDVI trend, or slope of time, is significant in both the univariate
model, and after accounting for climate variability in the multivariate models. Outcome B is a trend that
can be interpreted as driven by climate change and occurs where the NDVI shows a significant trend in
the univariate model but, this trend is not significant after accounting for variability in climate, within
the multivariate model. Outcome C, where the slope of time is not significant in univariate model, but it
is significant in multivariate model, is interpreted as driven by local change. A significant trend in the
multivariate analysis would suggest that the change is presumably due to direct factors. However, the
trend in the univariate analysis might not be significant because the variation in NDVI associated
indirectly with variation in climatic factors was masked by the influence of direct factors on NDVI.
Finally, Outcome D where the slope of time is not significant in either the univariate or multivariate
model, is interpreted as no change.
10
-------
EPA 600/R-25/337 I April 2025
Our Approach
Univariate Autoregression
NDVl, = 6 -e * Time-/.i,
Multivariate Autoregression
JVDVIt = fio + fiiRFt + faTmint + maxt + fi4VDPmaxt +
Comparison
Univariate time trend (0J vs. multiple time trend (p5)
Multivariate Autoregression
Time coefficient ((35) significance
S (p<0.05) NS (p>0.05)
Trend after No trend after
accounting f or c I i mate accountingforclimate
U"> Q
9 z
o c
to
rjj
>
o 2
o -
Al T3
Q. c
•— a'
(-
Z O
A
Local factors
driving NDVl trend
B
Climate driving
NDVl trend
C
Local factors driving
NDVl trend
D
No trend in NDVl
Figure 2-2. Graphical abstract for our approach in comparing NDVl trend between univariate and multivariate regressions.
The four outcomes are: Outcome A (local change): slope of time is significant from both univariate and multivariate models.
Outcome B (climate signal): slope of time is significant in univariate model but not in multivariate model. Outcome C (local
change): slope of time is not significant in univariate model, but it is significant in multivariate model. Outcome D (no
change): slope of time is not significant in both univariate and multivariate models.
11
-------
EPA 600/R-25/337 I April 2025
3 Analyses for Climatic Factors over 31 Years (1989-2019) in
Tillamook Basin
3.1 Description
Overall summary statistics were first calculated for each 30 m x 30 m pixel in the basin for all climatic
factors. Several analyses for precipitation were conducted on a monthly and seasonal basis, and for wet
(October through April) and dry (May through September) seasons. Trend analyses applying univariate
autoregression were conducted for each of the climatic parameters (e.g., precipitation versus time). The
trend direction for each significant pixel was then spatially mapped to identify any geographic patterns
of significance and trend direction.
The statistical summary of the climatic factors is given in Table 3-1. Ambient temperatures ranged
between -7.0°C (20.0°F) for minimum temperature and 31°C (87°F) for maximum temperature.
Statistical values for dewpoint temperature were close to that of minimum temperature. The average
monthly precipitation over 31 years ranged from 0.1 to 1,746.7 mm/month (Table 3-1). Potential
Evapotranspiration ranged from 0 to 142.6 mm/month (Table 3-1), which is much less than that of
precipitation, indicating high soil moisture content during the rainy season. The moisture index varied
between dry (-1) to moist condition (+1). Soil wetness can be further explored by the spatial distribution
of moisture index (see Section 3.3). The apparent temperature (AT) index for heat waves related to
human health.
Table 3-1. Monthly, overall, and seasonal summary statistics (minimum, mean, maximum, coefficient of variation) for the
climatic factors used in the study over 31 years for the Tillamook basin. Note: hPa = hectopascal.
Climate Factor - monthly
Minimum
Mean
Maximum
CV(%)
Minimum Temperature (°C)
-6.98
5.28
15.26
74
Maximum Temperature (°C)
-0.85
14.40
30.51
44
Average Temperature (°C)
-3.87
9.82
22.17
51
Dew Point Temperature (°C)
-6.55
5.78
14.58
68
Maximum Vapor Pressure Deficit (hPa)
0.01
7.62
32.59
69
Minimum Vapor Pressure Deficit (hPa)
0.00
0.48
4.95
118
Precipitation (mm)
0.12
252.34
1,746.72
89
Potential Evapotranspiration (mm)
0.00
43.74
142.57
51
Moisture Index (unitless)
-1.00
0.54
1.00
100
Apparent Temperature (°C)
-5.71
7.33
20.76
74
Precipitation
Overall
0.1
252.3
1,746.7
89.2
Wet Season (October-April)
0.4
78.2
480.4
88.9
Dry Season (May-September)
0.1
376.7
1,746.7
57.2
Potential Evapotranspiration (mm/month)
0.0
43.7
142.6
51.2
Moisture Index (unitless)
-1.0
0.539
1.000
100.4
Precipitation (Total mm/season)
Wet Season
758.8
2,653.9
5,837.2
29.7
Dry Season
78.4
391.2
1,124.4
39.8
12
-------
EPA 600/R-25/337 I April 2025
The spatial distribution of average monthly precipitation (Figure 3-1A) was higher in the north and south
part of the basin and variability (coefficient of variation; CV) of monthly precipitation was higher in the
east than in the west part of the basin (Figure 3B). The spatial distribution pattern of average
precipitation values within the basin was the same in both dry and wet seasons (Figures 3-1C and 3-1E),
where the highest precipitation occurs in the center of the basin and decreases toward western and
eastern bounds. The range of precipitation values for the wet season is 1,645 to 3,956 mm/wet season
(Figure 3-1C). The range of precipitation values for the dry season is 203 to 571 mm/dry season (Figure
3E). Wet season precipitation is about 6.2 times higher than dry season (Figure 3-1C versus Figure 3-1E).
Variability (CV) of wet season precipitation is less than 30%, with CV is higher at the eastern side, then
decreases gradually toward the western side. Variability (CV) of seasonal precipitation in Figure 3-1F,
shows that variability is less than 50% for dry precipitation; CV is higher at the eastern side, then
decreases gradually toward the western side of the basin. The CV of dry season precipitation range is 33-
43, which is higher than that of wet season precipitation (20-26) (Figure 3-1D). Variability for the dry
season is lower in the middle basin; highest variability is found in the eastern part and around and
within the estuary.
The spatial distribution of monthly precipitation trend values showed that the highest decreases in
monthly precipitation were in the middle part of the basin (Figure 3-2A) with the highest negative rate
of -0.251 mm/month and the lowest negative rate of -0.052 mm/month are on the eastern and western
bounds of the basin. The inset in Figure 3-2B is monthly rainfall from the South Fork
(GHCND:USR0000OSFK) for ground truthing, where the slope direction and its p values match our
analyses results. While the monthly average precipitation ranged from 0.4 to 780.4 for the wet season,
the range was between 0.1 and 1,746.7 mm/month for the dry season (Table 2-1) reflecting presence of
an extreme individual storm events in the dry season.
Trend values for both wet and dry seasons were the lowest in the western and eastern bounds of the
basin (Figures 3-2 C & D) encompassing estuaries within the basin. Decreasing wet precipitation trends
are higher than that of dry precipitation; trend values for the wet season ranged from -5.2 to -31.8
mm/season compared to that of dry season (-0.4 to -5 mm/ season).
The index of monthly moisture (Ml) denotes the relative wetness of the environment and ranged from
very dry (-1) to very wet (+1). Ml values are positive when precipitation (PPT) > PET and become
negative when PPT < PET. Summary statistics for monthly Ml are shown in Table 2-1. The relationship
between precipitation and PET is influenced by landcover and climatic factors. Removal of vegetation via
clearcut, for example, increased evaporation and more soil moisture loss. While precipitation decreased
between years (Figure 3-3) in a tree removal spot area (Figure 3-4), PET remained high. Water loss
following old growth tree removal increased because new vegetation growth consumes more water
than old growth. Higher evapotranspiration from clearcut areas in forested watersheds decreases the
amount of water yield to streams (Segura et al., 2020; French et al., 2023) in the northwest USA.
Drinking water shortage was reported in many coastal communities in Oregon in the years before 2023
(French et al., 2023). In comparison to areas in which trees remained intact (Figure 3-5) for the same
years, PET is lower than that in areas were trees were removed (Figure 3-3).
13
-------
EPA 600/R-25/337 I April 2025
mm/Dry Season
Mini"*™1
| aw -
:M7 ¦
[
HmW-443
¦w-erc
Figure 3-1. Monthly precipitation within the Tillamook basin over 30 years (1989-2019; n=372 months). (A) average
precipitation (mm/month), (B) coefficient of variation [CV] of monthly precipitation. Average wet season total precipitation
(C) (mm/season) reveals that wet average precipitation values are higher in the middle of the basin and lower within the
western and eastern bounds. (D) coefficient of variation [CV] of wet season precipitation. Average dry season total
precipitation (E) (mm/season). (E) (coefficient of variation [CV] of seasonal precipitation) shows that variability <50% for dry
precipitation; CV is higher at the eastern side, then decreases gradually toward the western side of the basin.
14
-------
EPA 600/R-25/337 I April 2025
Te*w*ik SMtri Fp>K k!*iv RF Apr - Nov 2W#
Trend direction and Probabilty
| Decreased Significantly
I Decreased Not Significantly
Wet Trend Value mm/season Pry Trend Value mm/season
•31.8 - -26.5
|*28.4--21.2
-21.1 -.16.8
-15.7--10.5
10.4 - -5 2
12
I .1
0 $ 10 20 Kik>n-*wr«
> . ». t |. i ) <
Figure 3-2. The temporal trend of monthly precipitation in 30 m x 30 m pixels over 30 years (1989 - 2019). Map A (trend
values mm/month of monthly precipitation) shows that trend values decreased within the basin and ranged from -0.052 to -
0.251 mm/month. Map B (direction and probability of monthly precipitation) shows that monthly precipitation decreased
significantly (p < 0.05) in the middle of the basin. Rainfall station, South Fork (GHCND: USROOOOOSFK), is in Map B within the
area of non-significant precipitation trend. Slope direction and slope probability of precipitation from the rainfall station is
same as of that from Map B (non-significant grey are in the map trend from PRISM data). Seasonal precipitation trend values
within the Tillamook basin over 30 years (n=30 years). Map C (slope value of wet precipitation, mm/season) shows that the
magnitude of wet trends on the eastern and western part of the watershed are less than those in the middle of the
watershed in both seasons. Map D (slope value of dry season precipitation, mm/season) reveals that dry season precipitation
decreased with years but with lower rate than that of wet season, all are within the low class of wet season (< 5
mm/season). Prediction in legend biplot is predicted NDVI.
Trend value (mm/month)
¦-0 251 - -0.169
¦-0.168 - -0.118
El-0.117 • 0,052
a » *» m mo w
in«vj a rvrwJrt n*d
15
-------
EPA 600/R-25/337 I April 2025
month
pet PPT
Figure 3-3. Monthly precipitation (PPT) and potential evapotranspiration (PET) for 2006, 2007, and 2008. Trees were removed
after 2006. The location of this pixel is marked by the cyan symbol in Figure 6. Pixel size is 30- x 30-rn.
16
-------
EPA 600/R-25/337 I April 2025
|11, v/ater
a 21 .Developed mixed
22,Developed Lowlntensity
23, Developed, Mixed Intensity
H24. Developed, H igh Intensity'
H31, Barren
41, Deciduous Forest
42,E vergreen Forest
H)43, Mixed Forest
52, Shrub/Scrub
71 ,G rassland/H erbaceo us
81, Pasture/Hay
~82, Cultivated Craps
90, Woody Wetlands
95, Emergent Herbaceous Wetlands
0.0
1989 1992 1995 1998 2001 2004 2007 2010 2013 2016 2019
Time (Year)
Trend o NDVI Pred
Figure 3-4. 2008 National Land Cover Database (NLCD; https://www.mrlc.aov/data Of Tillamook basin; the cyan colored
point is within a clearcut area. Landscape type in areas where the marked point (Point 25) is located was evergreen in 2001 -
2006> grass (71) in 2008-2016, then shrub and scrub (52) in 2019. Biplot of NDVI vs time (inset) shows a dip around 2007. The
overall decrease in NDVI over 31 years was significant.
17
-------
EPA 600/R-25/337 I April 2025
2006
2007
1000
900
£
C
800
O
E
700
1
£
600
c
500
o
400
'a.
300
u
4J
200
EL
100
0
' ' ¦- a
' fl ' v
100
80 ¦=¦
60
40
—i 1 1 1 i i 1 1 1 r
1 2 3 4 S 6 7 8 9 10 11 12
month
- PET - PPT
1000
.
900
£
C
800
O
S
700
E
E.
600
c
500
o
400
'a.
300
o
oj
200
a.
100
0
Hi
OL
¦C
0
E
1
E.
c
o
IS
20
2008
1000
900
800
700
600
500
400
300
200
100
0
• v _
.'x
f
/
-V-
100
80 -p
60
40
LiJ
a.
20
1 2 3
4 5 6 7 8 9 10 11 12
month
- PET - m- PPT
100
80 £¦
c
0
1
60
•'W
lu
40 Q-
20
1 2 3 4 5 6 7 8 9 10 11 12
month
- PET - m- PPT
Figure 3-5. Monthly precipitation (PPT) and potential evapatranspiration (PET) for 2006, 2007, and 2008. Trees were intact in
2006, 2007, and 2009. The location of this pixel is marked on the Google Earth map in Figure App-4. Pixel size is 30- x 30-m.
3.2 Seasonal Anomalies
The aim of using the anomalies is to determine precipitation deviation over the 31-year time period for
both wet and dry seasons. The dry anomaly was calculated by subtracting the per-month average dry
season precipitation from the monthly dry season precipitation, and the wet anomaly was calculated by
subtracting the per-month average wet season precipitation from the monthly wet season precipitation
(Figure 3-6 A & B). Generally, precipitation in the wet season decreased but increased for the dry season
over the 31 years (Figure 3-6 C & D), Figure 3-6 C & D shows trends only for positive seasonal anomalies
because negative values did not show temporal trends. As shown in Figure 3-6, anomaly dry
precipitation increased with years whereas the trend decreased for wet season. Precipitation for
seasonal anomalies over the 31-years was also mapped to visualize spatial distribution.
The number of years with seasonal anomalies was higher in the center than in western and eastern
bounds of the basin (Figures 3-7 A & B). The eastern and western part of the basin are dominated by the
years with negative anomalies in Figures 3-7 C & D. From the four figures, one can summarize that the
central 75% of the basin experienced more precipitation anomalies over the 30-yr period than did the
coastal plain or the eastern edge of Tillamook basin.
18
-------
EPA 600/R-25/337 I April 2025
CO
£
£
1500
1000
500
-500
-1000
TF
un
1500
1000
500
-500
-1000
B
MJ u yj MJ-<- - Lr Lr
1989 1993 1997 2001 2005 2009 2013 2017
Year
1989 1993 1997 2001 2005 2009 2013 2017
Year
Wet season anomly > 0
10
c
o
to
a!
5
1600
1400
trt
.c
1200
*->
c
o
£
1000
rv
800
E
600
400
200
0
1980 1990 2000 2010 2020
Figure 3-6. Seasonal precipitation (mm/season) over 30-years for the selected pixel. The horizontal line at 0 in both figures
represents the seasonal precipitation 31-year average. (A) wet season anomaly precipitation and (B) dry seasonal
precipitation anomaly. Seasonal precipitation anomaly > 0 values for wet and dry season. (C) Wet season precipitation
anomaly values and its trend (y = -7.2745x + 15167, R2 = 0.0194). (D) Dry season precipitation anomaly values and its trend (y
- 5.7265X -11335, R2 = 0.19).
19
-------
EPA 600/R-25/337 I April 2025
Annual Wet Anomaly > avg
B
min=1
2-7
8-12
13-16
17-19
20-28
max=29
Annual Dry Anomaly > avg
B
min=2
3-9
10 -12
13-15
16-18
19-25
max=26
Annual Wet Anomaly < avg
I
min=1
2-10
11 -13
14-17
18-22
23-29
max=30
Annual Dry Anomaly < avg
in
A
5 10 20 Kilometers
i i I i i i I
min=5
6-12
13 -15
B 16-18
19-21
22-28
max=29
Figure 3-7. Number of years with precipitation anomalies. (A) Wet season number of years > 30-year average (2,653.9
mm/season). (B) Dry season number of years > 30-year average (391.1 mm/season). (C) Wet season number of years < 30-
year average (2,653.9 mm/season), and (D) Dry season number of years < 30-year average (391.1 mm/season).
20
-------
EPA 600/R-25/337 I April 2025
Maps of seasonal anomaly trend values in Figure 3-8 show that the wet anomalies trend decreased
within most of the basin (Figure 3-8A), the dry anomalies trend increased within most of the basin
especially in the western bound (Figure 3-8B) of the basin. Estuaries experienced an increase in dry
anomalies (Figure 3-8B).
WetAnomly DryAnomly
¦-18--15 [ZI1-3
~-14--12 [ZI4-14
~-11 -0
~1 -14
Figure 3-8. Distribution of wet and dry seasonal anomaly trends. (A) Wet anomaly > 30-year wet season average. (B) Dry
anomaly > 30-year dry season average. The darker shaded areas denote higher trend values of seasonal precipitation
(amount) anomalies. Estuaries experiences more precipitation within the dry season but experience no change from the wet
season. Blank areas indicate no trend evaluation.
3.3 Derived Climatic Factors: Apparent Temperature Index, Potential
Evapotranspiration, and Moisture Index
The apparent temperature (AT) index for heat waves is organized into four groups based on human
health impacts (Smoyer-Tomic and Rainham, 2001). The lowest group ranges from 26.7°C to 32.2°C and
long exposure to this AT range will result in heat fatigue. The highest apparent temperature for the
Tillamook basin was 20.76°C which is lower than the lower limit of the lowest AT group suggesting that
there were no significant health impacts from heat waves in the basin during the period of record. The
range of AT within the temporal range of this study did not change from 1989-2013, and then decreased
in western Oregon (Nash and Christensen, 2017). Since the maximum AT value is less than the lowest
limit proposed by Smoyer-Tomic and Rainham (2001) for "fatigue possible with prolonged exposure
and/or physical activity (26.7°C < AT < 32.2°C)/' no further analyses were conducted for AT.
Annual average values for PET, precipitation, and Ml over 31 years are presented in Figure 3-9. The PET
ranged from 0 mm/month to 143 mm/month with an annual average of 44 mm/month. Annual average
of Ml ranged from 0.37 to 0.76, which falls in the wet to very wet scales per Feddema's (2005) wetness
21
-------
EPA 600/R-25/337 I April 2025
classification. For one pixel in a clear-cut, the relationship between PET, Ml and precipitation monthly
averages over 31-year is presented in Figure 3-3. The location of this pixel is shown in Figure 3-4. High
Ml results when precipitation is more than PET and occurred in January-April and September-December.
Drier moisture conditions (low Ml) result when PET exceeds precipitation, which occurred in May-
October. Correspondingly, Ml became positive, reflecting moist conditions and became negative, in dry
conditions. The spatial distribution of Ml in the basin corresponds with the annual average precipitation
pattern in Figure 3-1A.
Figure 3-9. PET annual average (A), precipitation annual average (B), and annual average of moisture index (C). The darker
shaded areas indicate higher average values. Higher annual PET are coinciding with that of less precipitation. Moisture index
values ranged from 0.37 to 0.76 which falls in the wet to very wet scales per Feddema (2005) wetness classification. The
distribution of moisture index in the basin corresponds with that of precipitation pattern. Pixel size is 30 m x 30 m.
22
-------
EPA 600/R-25/337 I April 2025
The monthly Ml average value is 0.54 with range values of -1 to +1 (Table 2-1), and the average values in
Figure 3-9C are all positive within the basin. The basin experiences negative Ml, yet these extreme
values disappear from maps of average values. Hence, the analyses were run per month for these three
derived climatic factors but at a larger scale (800 m x 800 m pixels) (Figure 3-10) for better presentation
of pattern and ranges of precipitation (Figure 3-10 left panels), PET (Figure 3-10 middle panels), and Ml
(Figure 3-10 right panels). Spatial distribution of precipitation from the 800 m x 800 m pixel scale and
within the basin is similar to that of 30 m x 30 m pixel scale (Figure 3-10, right panels). Precipitation
within the central part of the basin is higher than that in eastern and western bounds, but the area with
high precipitation shrinks in July and August (Figure 3-10 left panels). The middle panels in Figure 3-10
are PET, which was highest in western bound of the basin with similar spatial distribution in January and
December. Estuaries are within the high PET areas. Following January, the PET pattern shifts to the east
and becomes highest in eastern bound from February to September. The right panel in Figure 3-10 is the
Ml, where Ml indicates wetness index as moist to very wet for January through May (Feddema, 2005).
The Tillamook basin experiences dry conditions where low and negative Ml appeared in the eastern
bound of the basin in June and extended to the whole basin on July, August, and September. Moisture
Index in August was driest.
23
-------
EPA 600/R-25/337 I April 2025
January
- Prep PET Ml
.fa* & > n§
iw #"TX 48k
B299-356 I 112-15
357 -427 16-17 II0.93 - 0.95
428 - 458 18 -20 ^¦0.96-0.98
459 - 488
489 - 653
February ^
TS- Kir
1229 280
^^J281 - 368 | 113-18
369 -414 [ [19 - 20 | |0.87 - 0.91
^M415 - 461 21 - 23 ^^0.92 - 0.96
¦¦462 - 580
March jj...
J!*f\ *mn Aik
*M* 4BT
BS:2^r "" ^ir
^343 -377 ^^25-29
^¦378 -406 I 30 - 32 0.81 -0.86
^5407 - 499 ^H33 - 36 ^HO^ - 0.92
"®>> m- mm-
B'" *R8ty
197 -228 ilii
229 -242 I |36 - 41 ~ I ,ne,
243 " 254 I 142 - 44
255 - 342 ^^45 - 49 ^B°-75 - 0.86
-------
EPA 600/R-25/337 I April 2025
June
B-0.39
-0.32
0.01 -
0.34-
- -0.33,
- 0, dry
0.33, moist
0.66, wet
-------
EPA 600/R-25/337 I April 2025
Prep
I]52 - 79
L_|80 - 94
95-102
^¦103-109
¦¦110 -147
-0.33 - 0, dry
]0 - 0.33,moist
10.69 - 0.88, very wet
November
0.9 - 0.97, very wet
October
334 - 395
y
9-13
14 -15
16-18
|0.94 - 0.98, very wet
396 - 491
¦492 - 541
^M542 - 589
^¦590 -735
Figure 3-10. Spatial distribution of average (n=31 years) precipitation (PPT, mm/month), potential evapotranspiration (PET,
mm/month), and moisture index (Ml, unitless) per month. Scale of Ml (Feddema, 2005; 1.00 =>MI >= 0.66 very wet, 0.33 <=
Ml < 0.66 very wet, 0 <= Ml < 0.33 moist, -0.33 <= Ml <0 dry, -0.66 <= Ml < -0.33 semi-arid, -1.00 <=MI < -0.66 arid). Pixel size
is 800 x 800 m.
26
-------
EPA 600/R-25/337 I April 2025
4 Trend Analyses for NDVI and related Climatic Factors over 31 years
(1989-2019) in Tillamook Basin
4.1 Temporal Trend in NDVI and Climatic Factors over 31 years
From the univariate autoregression model, the NDVI changed significantly in 56% of the Tillamook basin
over the 31-year period (1989-2019) in the univariate analyses (Table 4-1), The direction of the trend
was predominantly positive with 75% of the overall significant changes observed in areas with an
increase in NDVI (Table 4-1). NDVI trend values (Figure 4-1A) and its significant values (Figure 4-1B) cross
the basin shows this dominant increase in NDVI is mostly concentrated on the western side of the basin
closer to the Pacific coast (Figure 4-1). NDVI trend decreased significantly mostly where clustered in
eastern part of the basin.
Table 4-1. Percent of total areas with significant increasing or decreasing monthly trend in NDVI and climatic factors from
univariate autoregression model from 1989 to 2019 in the Tillamook basin, Oregon. The last two rows are Percent of total
areas with significant increasing or decreasing trend in NDVI and from univariate and multivariate autoregression models
from 1989 to 2019 in the Tillamook basin, Oregon.
Variable
Increase(%)
Decrease (%)
Total(%)
NDVI
42.13
14.28
56.41
Maximum Temperature
3.51
5.14
8.65
Minimum Temperature
95.59
-
95.59
Average Temperature
49.72
-
49.72
Max Vapor Pressure Deficit
44.36
10.20
54.56
Min Vapor Pressure Deficit
94.34
-
94.34
Average Dew Point
0.11
23.72
23.83
Precipitation
-
31.96
31.96
Univariate
42.1
14.3
56.4
Multivariate
28.0
4.7
32.7
Figure 4-1. Univariate autoregression model results: (A) trend values range between -0.0024 to + 0.0022, (B) significant
temporal trend for the monthly NDVI from 1989 through 2019 in the Tillamook basin, Oregon. Green pixels indicate
significant (p < 0.05) increase in greenness (i.e., NDVI); red pixels indicate significant decrease in greenness.
27
-------
EPA 600/R-25/337 I April 2025
The proportion of the Tillamook basin with significant changes in climatic parameters ranged from
approximately 9% (monthly maximum temperature) to approximately 96% (monthly minimum
temperature) (Table 4-1), There was no portion of the Tillamook basin that had statistically significant
decreases in monthly minimum temperature, monthly average temperature, or monthly minimum
vapor pressure deficit (Table 4-1), Figure 4-2 shows the spatial patterns of change for climatic factors
used in the multivariate regression model; patterns differed among the four climatic factors. There were
no areas of the Tillamook basin that had statistically significant increases in monthly precipitation over
the 31 years (Table 4-1; Figure 4-2A). The decreasing trend in precipitation agrees with that of OCCRI
(2013) where precipitation declined by 0.18 inches per decade or 2 inches overall from 1901 - 2011.
Spatial distribution of dewpoint and mean temperatures, and minimum vapor pressure deficit (Figure
App-2) are not used in the multivariate model but are used to explain changes.
Figure 4-2. From univariate autoregression model: pixels (30 m x 30 m) with significant temporal trend (1989-2019) for
monthly: (A) precipitation (red indicates significant [p < 0.05] decrease [drier]); (B) maximum temperature (blue indicates
significant decrease [cooler], red indicates significant increase [warmer]); (C) minimum temperature (red indicates significant
increase [warmer night temperature]); and (D) maximum vapor pressure deficit (red indicates significant increase and blue
indicates significant decrease). In all maps, grey indicates no significant change.
28
-------
EPA 600/R-25/337 I April 2025
The areas of significant decrease in precipitation clustered in the middle of the basin (Figure 4-2A).
While the minimum temperatures increased significantly throughout the basin (Figure 4-2C), the
significant increase in maximum temperature is concentrated into two small areas by the eastern border
of the Wilson and Trask watersheds (Figure 4-2B). There was a significant decrease in maximum
temperature clustered within the Tillamook and Trask watersheds, and smaller clusters of decreasing
maximum temperatures scattered within the Wilson and Kilchis watersheds (Figure 4-2B). Significant
increases and decreases in maximum vapor pressure deficit exhibited a geographic dichotomy with
increases in the eastern portion of the Tillamook basin and decreases in areas closer to the Pacific coast
(Figure 4-2D).
4.2 Changes of NDVI including Climatic Factors
When climatic factors are included in the multivariate model, the extent of the area with significant
change over the 31-year period of record decreased from 56% to 33% (Table 4-1, Figure 4-3). The
intensity of green pixels from multivariate is less than that for univariate in Figure 15. The NDVI
increased significantly for 28% of the pixels in the multivariate autoregression (Table 4-1) whereases it
increased significantly in 42% of the pixels in the univariate autoregression (Table 4-1).
Figure 4-3. Significant temporal trend for the monthly NDVI from 1989 through 2019 in the Tillamook basin, Oregon. (A)
univariate autoregression model and (B) multivariate autoregression model results: pixels (30 m x 30 m) with green pixels
indicate significant (p < 0.05) increase in greenness (i.e., NDVI); red pixels indicate significant decrease in greenness.
4.3 Climatic Factors Contribution
Coefficients of time from the multivariate and univariate regression models were summarized in a
contingency table (Table 4-2). The spatial distribution of the four outcomes linking changing in NDVI to
that of climatic or direct factors is presented in Figure 4-4 and Table 4-2. Outcome D (no significant
trend/change in the NDVI in both analyses) was the most frequent (42%). The NDVI significance
apparently resulted from factors other than just climate, such as direct factors like wildfire, because the
trend was significant regardless of whether climatic factors were included in the analysis. A total of 33%
(Outcome A: significance in both analyses [31.3%] and Outcome C: significant in the multivariate but not
the univariate analysis [1.3%]) of all the pixels were significant for the NDVI trend in the multivariate
29
-------
EPA 600/R-25/337 I April 2025
analysis, reflecting the influence of the direct, local factors on NDVI rather than the selected climatic
factors. Nearly 25% of the pixels showed climatic factors as the cause of the significant NDVI trend
(Outcome B: the trend was significant in the univariate analysis, but not in the multivariate analysis due
to inclusion of climatic factors in the multivariate model).
Table 4-2. Frequency of outcomes in comparison of NDVI trends between univariate and multivariate analyses for 3,291,543
30 m x 30 m pixels in the Tillamook basin, Oregon. Outcomes A and C (bold) denote changes due to local factors; Outcome B
(italics) denotes changes due to climatic factors, and Outcome D denotes no change in greenness.
Significant-
Multivariate
Non-Significant-
Multivariate
Total-
Multivariate
Significant-
Univariate
Outcome A
31.39
Outcome B
25.00
56.39
Non-Significant -
Univariate
Outcome C
1.37
Outcome D
42.23
43.60
Total -
Univariate
32.76
67.23
100.00
Figure 4-4. The spatial distribution (left) of the four outcomes (right) from comparison of the NDVI trend significance between
univariate and multivariate autoregression for each pixel. Definition of the four outcomes is in the contingence table (right).
Both Outcomes A & C reflect changes due local factors, Outcome B reflects a climate factor, and Outcome D reflects no
change.
30
-------
EPA 600/R-25/337 I April 2025
5 Cause of Greenness Change
Four types of analytical applications were explored further to examine dramatic changes in NDVI with
sufficient information available (information from previous in-depth studies and Google Earth imagery
over the 31-year period) to infer the likely cause for the change. Two examples include significant NDVI
changes related to the four outcomes: tree mortality in the whole basin; and south flow corridor (SFC)
restoration project in Tillamook Bay, Oregon.
5.1 Tree Mortality, Insect Infestation
Tree mortality data were obtained from the US Forest Service survey (USDA Forest Service, 2016) where
yearly polygons of infested areas (data collected early July through September) were mapped from 1989
to 2019 in the Tillamook basin. The mortality variable 'DAM1' was used to represent the number of dead
trees per acre. Tree mortality polygons were rasterized to 30 m x 30 m pixels and synchronized with that
of NDVI. Pixels with yearly maximum NDVI and tree mortality were extracted to examine the association
between NDVI and tree mortality.
The association of insects with hosts is dependent on climatic factors and varies with the geographical
location (Colorado State Forest Service, 2014; Liebhold and Bentz, 2017) where insect infestation
increased due to reduced resistance in drought-stressed trees. Forested areas across the entire basin
suffered tree mortality and defoliation from pest infestations and diseases from 1991 to 2019 (Figure 5-
1). The cumulative tree mortality from 1991 to 2019 (Figure 5-1D) was distributed throughout the basin,
with more clusters located in the northern part.
31
-------
EPA 600/R-25/337 I April 2025
] VJ
1 \ v-f'
x -**¦ C. ' (
YTT*"~K "v c
B
£*J i" - \ -~v
ift
' 1
X > ^-v
.--'• i
J "V • :-' ,^-V
. r~> "/
JL^J
gp*" Sv^. ' '--V
C
, x .. I r.?s' ^
/ '-K$m
K ¦ r?.V s?
r.K. -
/ ** V - I
l :¦ ^-+w; .•» v.r
* -Tt: /:V&'4£(
3s . ,v.
/ V.M
Figure 5-1. New incidence of tree mortality for (A) 1991, (B) 2006, (C) 2013, (D) 1990-2019 cumulative tree mortality. Map
with color polygons over the years in Figure App-4.
Defoliation arid die-off trees have a significant, negative effect on greenness (Figure 5-2). The significant
decrease in NDVI is attributable to direct factors of insect infestations that can be influenced by changes
in climate factors. There was spatial overlap with tree mortality (Figures 5-1 and 4-2A) and a decrease in
monthly precipitation. Further, most of the tree mortality locations overlapped with the large area of
monthly nighttime temperature increases (Figures 5-1 and 4-2C). An increase of the monthly minimum
temperature reduces the likelihood of extreme low temperatures which would set back pest infestations
(Pasek and Schaupp, 1992). The monthly maximum vapor pressure trend that significantly increased was
mostly in eastern part of the basin and decreased toward the western part of the basin (Figure 4-2D).
Vapor pressure deficit is a major factor for tree mortality in the southwestern US (Williams et alv 2013)
and Pacific Northwest (Jarecke et a!., 2024). Vapor pressure increases with increasing temperature and
evapotranspiration rate increases as VPD increases. This chain enhances tree mortality. Significant
rising minimum temperature and VPD and the significant decreasing trend in precipitation in
Figure 4-2A coincide with tree mortality polygons in Figure 5-1.
32
-------
EPA 600/R-25/337 I April 2025
Figure 5-2. A significant decreasing trend (slope = -0.00108 p < O.OOOl) between maximum NDVI and tree mortality (DAM1)
from 1989 through 2019 in the Tillamook basin.
Eight sample locations across the broader tree mortality area, and across a range of times, were
explored further (Figure 5-3). In ArcMap and within the marked area in Figure 5-3, we recorded the
slope direction and probability of the climatic factors in Table 5-1. The table shows the significant
increasing (+S) or decreasing (-S) trends within the eight areas. Minimum and mean temperatures, and
minimum vapor pressure deficit increased significantly within all eight marked areas (except for Sample
7 where minimum temperature was not significant). In general, trends in maximum temperatures
increased but were not statistically significant. While dewpoint temperature and precipitation decreased
(not statistically significant), vapor pressure deficit (both minimum and maximum) generally increased.
NDVI vs time for a pixel in Figure 5-3 is within the tree mortality polygon where tree damaged was
observed in 2002, 2005, and 2013. Changes in NDVI trend in response to tree damage should be gradual
indicating the abrupt NDVI change around March 2009 is in response to factor other than insect
infestation. Google Earth historical images identified a clear cut that took place between 2005 and 2011.
Vegetation regrowth was observed in 2012, and vegetation thickness increased by 2020. The slope of
NDVI before March 2009 was not significant and NDVI began to increase significantly (p<0.0001) after
March 2009 with slope = 0.0041.
33
-------
EPA 600/R-25/337 I April 2025
3/2009
>
Q
0.8
0.6
0.4
0.2
O ° «r q00o 00 0*o ^
goo 0 9 o™ o °^m I
n » ° n n n o
° °
o O o o
& £
>;8
I V
I s°%
I
b °
100
200
Time (month)
300
Figure 5-3. Locations of tree mortality samples. Marked polygons include tree damage from one year to several. Tree
mortality observed in Sample 1:1996; Sample 2:1991,1993; Sample 3:1991,1992, 2006, 2016; Sample 4:2019; Sample 5:
1991,1992, 2004, 2006, 2015; Sample 6:1993; Sample 7:1991, 2006, 2013, 2015; Sample 8:2013. The inset biplot shows the
NDVI temporal behavior in a pixel (row=1692 col=210) within Sample 8, with abrupt change around 2008-2009.
34
-------
EPA 600/R-25/337 I April 2025
Table 5-1. Univariate trend for climatic factors for Sample Locations (SL) within a tree damaged sample (pixel: 30 m x 30 m).
Tree damage was observed in sample (year): Sample 1 (1996); Sample 2 (1991,1993); Sample 3 (1991,1992, 2006, 2016);
Sample 4 (2019); Sample 5 (1991,1992, 2004, 2006, 2015); Sample 6 (1993); Sample 7 (1991, 2006, 2013, 2015); Sample 8
(2013). Vertical shaded columns denote that tree damaged was observed in more than one year.
Climatic Factor
SL #1
SL #2
SL #3
SL #4
SL #5
SL #6
SL #7
SL #8
Min-Temperature
+s*
+s
+s
+s
+s
+s
+s
+s
Max-
Temperature
NS
NS
NS
NS
NS
NS
NS
NS
Mean
Temperature
+s
+s
+s
+s
+s,
+s
NS
+s
Dewpoint
Temperature
-s
-s
-s
NS
NS
NS
NS
NS
Vapor Pressure
Deficit -
minimum
+s
+s
+s
+s
+s
+s
+s
+s
Vapor Pressure
Deficit -
maximum
+s
+s
+s
+s
+s
-NS
-s
-s
Precipitation
NS
-s
-s
NS
NS
NS
NS
NS
*+S trend is positive significant, -S trend is negative significant, NS trend is not significant.
Note: Information in the above table was extracted manually from overlaying climatic factors and tree
mortality layers in ArcMap.
5.2 Changes in Greenness in Restoration Areas
We utilized a restoration polygon in the SFC, Tillamook County, Oregon (Janousek et al., 2021; Figure 5-4
bottom) to follow changes between pre- and post-restoration in specific locations to explain our
analyses outcomes that relate to local factors. Landcover/landuse in the restoration area is diverse and
includes forested, crop, and grazing areas. Elevation is higher at northern side than the south, the grazed
zone is at a higher elevation than the cropped zone. Prior to restoration, this tidal wetland was diked for
many decades and was used for agriculture and grazing. But several years before restoration, agriculture
was abandoned. Restoration started by removing dikes, and water moved through the corridor along
the slope gradient filling and widening channels.
We posted the restoration polygons provided by Laura Brophy (Janousek et al., 2021) on our 4 outcomes
map and chose a sample within Outcome D and another sample within Outcome A (Figure 20) to
compare vegetation changes with that of NDVI. Although the vegetation samples in pre-restoration
(2013-2015) and early post-restoration (2017-2020) fall in relatively narrow window within our study
periods (1989 - 2019), they still can be utilized to ground truth our results.
Outcome D sample is within non-tidal forested wetlands (green) polygons in Figure 5-4 and Figure App-3
where the dominant vegetation was Sitka spruce (Picea sitchensis). Sitka spruce density did not change
from pre- to post- dike removal due to presence of fresh water which was not influence by inundation
(Janousek et al., 2021). The NDVI temporal trend over 31 years was not significant (Figure App-3) and
did not change upon including climatic factors in the multiple regression (p>0.06), therefore this sample
represents an Outcome D.
35
-------
EPA 600/R-25/337 I April 2025
The sample within Outcome A is located in the south zone of the south flow corridor (Figure 5-4). NDVI
temporal trend (Figure App-5) decreased significantly in univariate regression (p<0,0003) and remained
significant in multivariate multiple regression (p<0.017); the trend was not influenced by including
climatic factors. The dominant vegetation in this location is non-native vegetation, creeping bentgrass
(green pixel denotes Outcome A, Figure 5-4) which existed in post- and pre-restoration. However, native
vegetation decreased from 2014 to 2018 (Janousek et al., 2021).
A
Figure 5-4. Top-right, location of the restoration polygon within the Tillamook Basin; top-left, green, orange, and grey
polygons denote native dominant vegetation, non-dominant vegetation, and bare ground, respectively, as described by
Janousek et al. (2021). Bottom, selected restoration polygons from Janousek et a I. (2021) overlay two of the four outcomes (A
& D). Outcome D sample is within the northern zone of the SFC and Outcome A is in the southern zone of the SFC. Elevation
gradient sorted as Outcome sample D is at higher elevation than sample A. The locations of the outcome samples on Google
Earth and the biplotofNDVI vs months are in Figures App-3 and App-4 for Outcome D sample. Biplot of NDVI vs months for
Outcome sample A is in Figure App-5.
Wilson
Outcome: D
Outcome
0 12 4 Kilometers
1 ' i i | i i i |
36
-------
EPA 600/R-25/337 I April 2025
5.3 Changes in Greenness due to Climatic Factors
We selected a sample (one raster) that is within Outcome B. This sample is in the Trask watershed and
within Outcome B where climatic factor influenced greenness (Figure 5-5). This raster was within mixed
forest between 1989 - 2019; the landcover type did not change over the study period at this location.
However, the greenness decreased significantly in univariate autoregression, and that change became
not significant when climatic factors were included in the multivariate regression. Tree condition and
mortality from insect infestation typically responds to climatic factors slowly over time. Hence changes
in greenness at sampled raster are likely attributable to climatic factors. East to the samples area (Figure
5-6), US Forest Service survey (USDA Forest Service, 2016) mapped tree mortality polygon that occurred
in 2008. NLCD showed changes in landcover type in the west side of this sample location in 2008, from
evergreen to herbaceous and shrub and scrub (Figure 5-6).
Figure 5-5. Outcomes related to climatic factors (top, green denotes local changes, dark red denotes changes related to
climatic factors).
37
-------
EPA 600/R-25/337 I April 2025
Developed mixed
Deciduous Forest
|Evergreen Forest
Mixed Forest
Shrub/Scrub
Grassland/Herbaceous
Pasture/Hay
Woody Wetlands
Figure 5-6. Historical NLCD for areas around sample raster (Cyan closed circle). This sample remains within a tree intact area
in 1994, 2003, 2005, 2011, 2016, 2019, and 2024, whereas tree removal to the west of the sample location took place
between 2005 and 2011.
38
-------
EPA 600/R-25/337 I April 2025
6 Major Findings and Conclusions
6.1 Major Findings
Leveraging temporal and spatial climatic factors data helped in determining changes in climatic factors
over 31 years. Precipitation decreased significantly (p < 0.05) over the 31 years in nearly one third of the
basin. For one pixel in this area, precipitation decreased from 391 mm to 342 mm between the
beginning and end of this study. Following precipitation, average dewpoint decreased in 24% of the
basin and max vapor pressure deficit (VPD) decreased in 10% of the basin. Night temperature and
minimum vapor deficit increased in more than 94% of the basin, whereas maximum vapor pressure
deficit and mean temperature increased in nearly 50% of the basin (Table 4-1).
Localized heterogeneity in monthly PET and Ml spatial distribution maps were observed (Figure 3-9).
Hence, the analyses for PET and Ml were done at a coarser spatial structure (800 m x 800 m) (Figure 3-
10) to obtain a clearer cluster pattern. The Ml ranged from moist to very wet in September through April
but reached the arid level in July and August. The analysis of 31-year monthly averages of climatic
factors analysis showed that the basin went through dryness in June-September and extreme dryness in
July and August when the moisture index was in the arid scale of Ml (Figure 3-10). The spatial
distribution of PET exhibited a gradual movement from the west to east, January through August, then
back to west.
Precipitation trend values for both wet and dry seasons were the lowest in the western and eastern
bounds of the basin, encompassing the estuaries within the basin. The decreasing wet precipitation
trend is higher than that of dry precipitation; trend values for the wet season ranged from -5.2 to -31.8
mm/season compared to that of dry season (-0.4 to -5 mm/ season).
Outcomes from univariate and multiple autoregression models indicated that 42% of the basin
experienced no significant changes in greenness directly or indirectly from 1989 to 2013. One third
(33%, Outcomes A + C) of the basin experienced changes due to direct or local factors and are mainly
clustered within the Tillamook and Trask watersheds. Several small clusters of pixels were affected by
climatic factors but not local factors, comprising 25% of the basin and mainly in the Miami and Trask
watersheds. The four outcome results (Figure 5-5) appear to be clustered in small patches reflecting the
finer scale of the NDVI data used, similar to the observation which was reported in Nash et al. (2017).
Using remote sensing data (e.g., NDVI) and PRISM climatic factors for long-term monitoring to
characterize changes in climatic, greenness, and correlation between them is a simple and relatively
inexpensive method. Our methodology detected broad-scale changes (slow and rapid) that related to
climatic or local influence. This kind of monitoring can provide decision-makers with early warning
signals for widespread general trends as well as a means to identify specific areas where land conditions
are degrading. However, to determine the cause of greenness change in an area requires knowledge
about the changes in landscape, climatic factors, and environmental condition. Our demonstration of
NDVI pattern changes in selected areas as related to certain causes, provide a way for users to
determine the cause. An increase in NDVI in agricultural areas often occurs in response to expansion of
cultivated areas (e.g., Tillamook watershed). This is the same in wetland restoration, where NDVI
increases as vegetation cover increases (e.g., southern flow corridor, Tillamook County, Oregon;
39
-------
EPA 600/R-25/337 I April 2025
Janousek et al., 2021). A decrease in NDVI occurs when cultivation in agriculture areas is abandoned
(Figures 3-1B; see also S4 in Nash et al., 2014) and in areas with insect infestations of trees (Colorado
State Forest Service, 2014; Figure 7 in Ash et al., 2017). Closer examination of the biplot of NDVI vs time
in any of the above areas, one can locate the approximate date when change took place. The trend of
NDVI before and after an event can be determined to describe direction and significance of change.
The comparison between two regression models, inclusions of climatic factors in the multiple regression
and their effect on the trend of time provided a way to distinguish between changes that directly or
indirectly related to changes in climatic factors over the 31 years. The result of model outcomes is
dependent mainly on the on the data used in the analyses. Our statistical approach allows for time
series data where it accounts for serial correlation. The model inspects 120 lags (120 months) for any
possible serial correlation.
Spatial and temporal presentation of data gave different representations for range values for climatic
factors from mapping overall moisture index (Ml) average values for the basin (Figure 3-9). For better
presentation of Ml, yearly averages (n=31 years) of precipitation (PPT), potential evapotranspiration
(PET), and Ml per month were mapped in Figure 3-10. Monthly maps for these three factors showed
how PET shifted from west to east in the basin and that the wetness index (Ml) was very wet for January
through May. The Tillamook basin experienced dry conditions where low and negative Ml appeared in
the eastern bound of the basin in June and extended to the whole basin in July, August, and September.
Ml in August was driest. Identification of ranges of Ml from negative to positive over space and time was
not possible when the overall average values were used (Figure 3-9).
6.2 Conclusions
This study presents a cost-effective method analogous to that in long-term ecological research to
monitor and characterize long term changes in greenness. Analyses of long-term data on climatic factors
and greenness can potentially detect broad-scale, slow changes, such as those caused by climate change
over decades, as well as more local and rapid changes such as those caused by fire, agriculture, land
clearing, and habitat restoration. Selecting specific areas, leveraging study results and expert knowledge,
we tracked change in greenness due to anthropogenic and natural stressors on different scales (from
minimum scale, 30 m x 30 m, to clusters of pixels).
Monitoring is an essential element in natural resources management. Conservation of natural resources
requires a well-established plan that includes a monitoring program to identify and validate outcomes.
These analyses can provide environmental decision-makers with early warning signals for widespread
general trends as well as to identify specific areas where land conditions are degrading or improving.
Subsequent to identifying areas classified by the four outcomes, one can evaluate the cause of changes
such as fire or fire recovery, clear cut or recovery from clear cut, change in agriculture extent or practice,
tree mortality due to insect infestation, and others.
Our intent in this report is to make the methodology, data, and analyses available for research efforts to
better understand the trends of NDVI and climatic factors and their interactions in hydrological, societal,
biological, and ecological systems within the Tillamook Basin. Increasing trend in temperature in
estuaries, for example, has a direct effect on natural resources, higher temperature enhances algal
blooms and lower oxygen levels.
40
-------
EPA 600/R-25/337 I April 2025
7 Quality Assurance
This study was conducted under the approved Quality Assurance project plan "L-PESD-0034581-QP-1-0,
Tillamook Normalized Difference Vegetation Index (NDVI)." Datasets used are described in the Materials
and Methods section.
We present a means of long-term monitoring using spatial and temporal greenness index (NDVI) and
climatic factors (PRISM) in conjunction with statistical models to determine ground condition. This kind
of study can be updated continuously by adding temporal data. This method has to be mutually inclusive
with field expert knowledge to make correct decisions and actions to improve or keep the status of land
condition. An additional point is the scale of the factors (or variables) used in the analyses. Previously 8
km x 8 km and 1 km x 1km pixels were used in many landscape models. This study used 30 m x 30 m
pixels for NDVI and downscaled climatic factors from 800 m x 800 m to 30 m x 30 m pixels.
The abrupt changes in temporal NDVI in a pixel caused by tree removal in clear cut areas are easily
spotted from our NDVI analyses and were validated from NLCD and historical images in Google Earth
map. Changes in NDVI related to local factors (Outcomes A and D) were also identified and validated
within restoration polygons where changes in vegetation were recorded before and after restoration.
41
-------
EPA 600/R-25/337 I April 2025
8 References
Abatzoglou, J.T., Williams, A.P. 2016. Impact of anthropogenic climate change on wildfire across western
US forests. Proceedings of the National Academy of Sciences, 113(42): 11770-11772.
Alashan S. 2020.Combination of modified Mann-Kendall method and Sen innovative trend analysis.
Engineering Reports, 2: el2131.
Anyamba, A., Tucker, C.J. 2005. Analysis of Sahelian vegetation dynamics using NOAA-AVHRR NDVI data
from 1981-2003. Journal of Arid Environments, 63: 596-614.
Bernal A.A., Kane J.M., Knapp E.E., Zald, H.S.J. 2023. Tree resistance to drought and bark beetle-
associated mortality following thinning and prescribed fire treatments, Forest Ecology and
Management, 530: 120758.
Bonfils, C., Santer, B., Fyfe, J., Marvel, K., Phillips, T., Zimmerman, S. 2020. Human influence on joint
changes temperature, rainfall, and continental aridity. Nature Climate Change. 10: 726-731. URL:
https://doi.org./10.1038/s41558-020-0821-l.
Chhabra, A., Dadhwal, V.K. 2004. Estimating terrestrial net primary productivity over India using satellite
data. Current Science, 86(2): 269-271.
Colorado State Forest Service. 2014. Colorado forest insect and disease update. A supplement to the
2014 report on the health of Colorado's forests. URL:
https://csfs.colostate.edu/media/sites/22/2015/03/Final-2014-lnsect-Disease-Update-
2March2015.pdf.
Deininger K., Ali, D.A., Kussul, N., Shelestov, A., Lemoine, G., Yailimova, H. 2022. Quantifying war-
induced crop losses in Ukraine in near real time to strengthen local and global food security.
Development Economics, Development Research Group, Policy Research Working Paper. 10123.
URL: https://documents.worldbank.org/en/publication/documents-
reports/documentdetail/099102207072236777/idu0de02933f001de04df5087c30538da5ba4b35.
Feddema, J.J., 2005. A revised Thornthwaite-type global climate classification. Physical Geography,
26(6): 442-466.
French, E., Edulbehram, U., Hughes, S., Arndt, M. 2023. Oregon Coast Range ecological conservation:
Mapping recent logging within drinking watersheds of Oregon's Coastal Range to support future
resource management policies. Technical Report August 11, 2023. URL:
https://ntrs.nasa.gov/citations/2023001219Q.
Janousek C., Bailey, S., van de Wetering, S., Brophy, L., Bridgham S., Schultz, M., Tice-Lewis, M. 2021.
Early post-restoration recovery of tidal wetland structure and function at the Southern Flow
Corridor project, Tillamook Bay, Oregon. Oregon State University, Tillamook Estuaries Partnership,
Confederated Tribes of Siletz Indians, Institute for Applied Ecology, and University of Oregon. URL:
https://appliedeco.org/wp-content/uploads/Janousek et al 2021 SFC monitoring rpt.pdf.
Jarecke K.M., Bladon K.D., Meinzer F.C., Wondzell S.M. 2024.lmpact of rainfall and vapor pressure deficit
on latewood growth and water stress in Douglas-fir in a Mediterranean climate. Forest Ecology and
Management, 551: 121529.
42
-------
EPA 600/R-25/337 I April 2025
Ji, L., Peters, A.J. 2003. Assessing vegetation response to drought in the northern Great Plains using
vegetation and drought indices. Remote Sensing of Environment, 87: 85-98.
Konings A.G., Piles, M., Das, N., Entekhabi, D. 2017. L-band vegetation optical depth and effective
scattering albedo estimation from SMAP. Remote Sensing of Environment, 198: 460-470.
Koopman M.E. 2018. Climate change preparedness strategy for Tillamook Estuaries Partnership. GEOS
Institute. URL: https://climatewise.org/wp-content/uploads/proiects/tep-climate-change-
report.pdf.
Land Degradation Assessment in Mediterranean Europe (LADAMER). 2003. Land degradation
assessment in Mediterranean Europe. EC DG Research. Contract No EVK2-CT-2002-00179.
Contribution to the GMES Report. URL: https://www.copernicus.eu/en/land-degradation-
assessment-mediterranean-europe.
Liebhold, A., Bentz, B. 2017. Insect disturbance and climate change. U.S. Department of Agriculture,
Forest Service, Climate Change Resource Center. URL: https://www.fs.usda.gov/ccrc/topics/insect-
disturbance-and-climate-change.
Lybbert, T.J., Aboudrare, A., Chaloud, D.J., Magnan, N., Nash, M.S. 2011. Booming markets for Moroccan
argan oil appear to benefit some rural households while threatening the endemic argan forest.
Proceedings of the National Academy of Sciences, 108: 13963-13968.
Miller, D.L., Alonzo, M., Meerdink, S.K., Allen, M.A., Tague, C.L., Roberts, D.A., McFadden, J.P. 2022.
Seasonal and interannual drought responses of vegetation in a California urbanized area measured
using complementary remote sensing indices. Journal of Photogrammetry and Remote Sensing, 183:
178-195.
Najafi, M.R., Zwiers, F.W., Gillett, N.P. 2016. Attribution of the spring snow cover extent declines in the
Northern Hemisphere, Eurasia and North America to anthropogenic influence. Climate Change,
136(3-4): 571-586.
Nash, M.S., Bradford, D.F., Wickham, J.D., Wade, T.G. 2014. Detecting change in landscape greenness
over large areas: An example for New Mexico, USA. Remote Sensing Environment, 150: 152-162.
Nash, M.S., Wickham, J., Christensen J., Wade, T. 2017. Changes in landscape greenness and climatic
factors over 25 years (1989-2013) in the USA. Remote Sensing, 9(3): 295.
Nash, M.S, Christensen, J.R. 2017. Description of changes in climatic indices in USA over 25 Years (1989 -
2013). EPA/600/S-17/256. URL: https://www.epa.gov/enviroatlas/description-changes-climatic-
indices-usa-over-25-vears-1989-2013.
OCCRI (Oregon Climate Change Research Institute). 2013. Climate Change in the Tillamook Bay
Watershed. Oregon Climate Change Research Institute. URL:
https://digitalcollections.librarv.oregon.gov/nodes/view/166834.
Pasek, J.E., Schaupp, W.C. Jr. 1992. Populations of Douglas fir beetle in green trees three years after the
Clover Mist Fire on the Clarks Fork Ranger District. Shoshone National Forest, Wyoming; U.S. Forest
Service Rocky Mountain Region, Renewable Resources Staff: Denver, CO, USA, 1992; pp. 13. URL:
https://www.biodiversitvlibrarv.Org/item/293191#page/5/mode/lup.
43
-------
EPA 600/R-25/337 I April 2025
Perry, A.G., Korenberg M.J., Hall, G.G., Moore, K.M. 2011. Modeling and syndromic surveillance for
estimating weather-induced heat-related illness. Journal of Environment and Public Health, 1-10.
Peters, A.J, Eve, M.D., Holt, E.H., Whitford, W.G. 1997. Analysis of desert plant community growth
patterns with high temporal resolution satellite spectra. Journal of Applied Ecology, 34: 418-432.
PRISM. 2021 Descriptions of PRISM spatial climate datasets for the conterminous United States. Last
revised November 2021. URL: https://prism.oregonstate.edu/documents/PRISM datasets.pdf.
Ruggiero, P., C.A. Brown, P. Komar, J. Allan, D. Reusser, Lee, H. 2010. Chapter 6. Impacts of climate
change on Oregon's coasts and estuaries, Chapter 6 in "Oregon Climate Change Assessment Report".
U.S. Environmental Protection Agency, Washington, DC, EPA/600/R-10/184/2011. URL:
https://pnwcirc.org/sites/pnwcirc.org/files/ocar2010.pdf.
Segura C., Bladon, K.D., Hatten, J.A. Jones, J.A., Hale, V.C., Ice, G.G. 2020. Long-term effects of forest
harvesting on summer low flow deficits in the Coast Range of Oregon. Journal of Hydrology, 585:
124749.
Smoyer-Tomic, K.E., Rainham, D.G.C. 2001. Beating the heat: development and evaluation of a Canadian
hot weather health-response plan. Environmental Health Perspectives, 109(12): 1241-1248.
SAS Institute Inc. 2016. SAS/ETS 14.2 User's Guide. Cary, NC: SAS Institute Inc. URL:
https://documentation.sas.eom/doc/en/etscdc/14.2/etsug/main contents.htm.
TBNEP (Tillamook Bay National Estuary Project). 2020. 2020 State of The Bays. Garibaldi, OR. URL:
https://media.iournoportfolio.eom/users/76420/uploads/8dblac98-940d-4467-a308-
b87ba85c6b06.pdf.
Thorne, J., Boynton, R., Flint, L., Flint, A., Le, T. 2012. Development and application of downscaled
hydroclimatic predictor variables for use in climate vulnerability and assessment studies. California
Energy Commission. CEC-500-2012-010. URL:
https://ca.water.usgs.gov/pubs/Flint BMC California.pdf.
Tucker, C.J., Anyamba, A. 2011. Historical Perspectives on AVHRR NDVI and Vegetation Drought
Monitoring, https://core.ac.uk/download/pdf/10561646.pdf.
Turvey, C.G., Mclaurin, M. 2012. Applicability of the Normalized Difference Vegetation Index (NDVI) in
index-based crop insurance design. Weather, Climate and Society, 271-284.
Twumasi, Y.A., Coleman, T.L., Manu, A., Merem, E.C., Osei, A. 2011. Relationship between climate
parameters and forest vegetation at and near Digya National Park, Ghana. British Journal of
Environment and Climate Change, 1: 201-215.
USDA Forest Service. 2016. Aerial Insect and Disease Survey GIS Data for Oregon and Washington 1947-
Present. URL: https://www.fs.usda.gov/detail/r6/forest-grasslandhealth/insects-
diseases/?cid=stelprd3791643.
Vanderhoof M.K, Lane, C.R., McManus, M.G., Alexander, L.C., Christensen, J.R. 2018. Wetlands inform
how climate extremes influence surface water expansion and contraction. Hydrology and Earth
System Sciences, 22, 1851-1873.
44
-------
EPA 600/R-25/337 I April 2025
Wang, J., Rich, P.M., Price,K.P. 2003. Temporal response of NDVI to precipitation and temperature in the
central Great Plains, USA. International Journal of Remote Sensing, 24: 2345-2364.
Weiss, J.L., Castro, C.L., Overpeck, T.J. 2009. Distinguishing pronounced droughts in the southwestern
united states: seasonal and effects of warmer temperatures. Journal of Climate, 22: 225918-5932.
Williams A.P., Allen, C.D., Macalady, A.K., Griffin, D., Woodhouse, C.A., Meko, D.M., Swetnam, T.W.,
Rauscher, S.A., Seager, R., Grissino-Mayer, H.D., Dean, J.S., Cook, E.R., Gangodagamage, C., Cai, M.,
McDowell, N.G. 2013. Temperature as a potent driver of regional forest drought stress and tree
mortality. Nature Climate Change, 3: 292-297.
Williams, P., E.R. Cook, J.E. Smerdon, B.I. Cook, J.T. Abatzoglou, K. Bolles, S.H. Baek, A.M. Badger, and B.
Livneh. 2020. Large contribution from anthropogenic warming to an emerging North American
megadrought. Science, 368: 6488.
45
-------
EPA 600/R-25/337 I April 2025
Appendix 1
2016
2019
2006
2008
2013
¦Water
~Developed mixed
]Developed Low Intensity
B Developed, Mixed Intensity
Developed, High Intensity
~Barren
Deciduous Forest
(¦Evergreen Forest
~Mixed Forest
~Shrub/Scrub
~Grassland/Herbaceous
~Pasture/Hay
Cultivated Crops
~Woody Wetlands
~Emergent Herbaceous Wetlands
Figure App-1. Landcover/landuse from 2001 to 2019 in Tillamook basin, Oregon.
2004
46
-------
EPA 600/R-25/337 I April 2025
The pattern of increasing vapor pressure deficit (Figure App-2C) is similar to that of minimum
temperature (Figure 4-2C). The spatial distribution of significant trends for average and dewpoint
temperatures are mostly in the eastern portion of the basin. Significant increasing trend of average
temperatures are at the ridge top between watersheds. This pattern is similar in dewpoint temperature
in these areas. Average temperature (°C), shown in Figure App-2A, increased significantly in 50% (Table
4-1) of the basin, dewpoint temperature (°C), shown in Figure App-2B, decreased in 24% of the basin,
and minimum vapor pressure (hPa) increased in 94% of the basin. The pattern of increasing vapor
pressure deficit (Figure App-2-C) is similar to that of minimum temperature (Figure 4-2C). The spatial
distribution of significant trends for average ambient temperature and dewpoint temperature are
mostly in eastern bound of the basin. Significant increasing trends of average temperatures are at the
ridge tops between watersheds. This pattern is similar in dewpoint temperature in less areas.
C
A
0 5 10 20 Kilometers
1 i i i I i I i I
Figure App-2. Monthly significant trend within the Tillamook basin over 30 years (1989-2019; n=372 months) for average
temperature (°C) (A), dewpoint temperature (°C) (B) and minimum vapor pressure deficit (hPa) (C). Red shaded areas denote
significant increase, shaded green areas denote significant decrease.
47
-------
EPA 600/R-25/337 I April 2025
1989 1995 2001 2007 2013 2019
Time (Year)
Figure App-3. Google Earth map showing location of sample for an Outcome D in 2005 (A) and 2011 (B). Trend of monthly
NDVI (C) for the sample in (A) and (B) decreased non significantly (p=0.057) and remained not significantly decreasing when
climatic factors were included (P=0.07).
48
-------
EPA 600/R-25/337 I April 2025
/ v-
. :;h
" -v %, M ^
W' i' -tf ¦ T,.Vr
K * - V * . . . v V'-r' /• vj»->
i * ¦( - - — .
Tillamook Bay i • . • • . ^ '«* ¦/"--<¦¦ " ^
* »>* yj" xK-**'
"
• \ •*»• ••' • • ;S<- " >
L ¦•?, J < .... - v u>^- r^-\
V-v-i'V-- ', *"
i
„ 'a. >'-T
«s " '¦ ' -v*-v--"a ¦'— ¦ F\ /\ ' ~ . V M V •' ^
" , • '¦¦ .y ;' ¦¦ .
4;' N X. 4 , Tr«K River • • ' *4 V'l ' ' >'
/ " "Xx ~ * ** ' i|% • •"* >f -. ¦", ;**'•>
V- 1 •
y. .• V tf-". ..- :;"V> • V .v./
,
-------
EPA 600/R-25/337 I April 2025
Slope=-0 0064 p= 0.0003
1 1 1 1 1
1989 1995 2001 2007 2013 2019
Time (Year)
Figure App-5. Google Earth map for Outcome A sample location and its NDVI vs months.
50
-------
EPA 600/R-25/337 I April 2025
Appendix 2
Landsat collection 1, tier-1 level time series data used in this study were hosted in the Google Earth
Engine (GEE) cloud-computing platform data catalog. The Earth Engine (EE) code editor environment
was used to generate 30-meter resolution NDVI time series monthly maximum value composites for
available satellite data scenes between the years of 1989 to 2019 that corresponded with the Tillamook
Watershed boundary extent. Given the lengthy proposed study timeline and the limited availability of
image scenes afforded by the Landsat satellite orbit 16-day location revisit cycle,
interpolation/smoothing of data was not incorporated for the initially calculated monthly NDVI
maximum value composites. Values -2, -3 were excluded from the data. Proc Expand in SAS 9.2 (SAS®,
Cary, North Carolina) with method "Join" was used to substitute for missing values. The join function
integrates linear interpolation lines of non-missing values over all NDVI values. Interpolated values are
then used to fill the missing.
Landsat 8 Operational Land Imager (OLI), Landsat 7 Enhanced Thematic Mapper Plus (ETM+) and
Landsat 4-5 Thematic Mapper (TM) calculated Top-of-Atmosphere (TOA) reflectance data products were
matched for corresponding study timelines. Landsat image tiles whose metadata collection date
timestamps corresponded with the study timeline were selected from within the GEE data catalog for
subsequent use in the study. As this study focused on general long-term trend mapping, TOA product
collections, albeit uncorrected for atmospheric effects, contain a greater number of total available
image scenes as compared to Surface Reflectance (SR) product derivatives whose atmospheric data
correction processing were impacted by missing auxiliary input data and/or necessary thermal data.
There are a few established caveats pertaining to the direct use of TOA reflectance products in
timeseries analyses. There are cross-platform spectral response differences, particularly between
Landsat 7 ETM+ and Landsat 8 OLI sensors, attributed to the inherent differences in radiometric
resolution between sensors. These inherent differences result in a higher frequency of saturated pixel
values across spectral bands in Landsat 7 ETM+ data products. Due to the radiometric calibration
differences, Landsat 7 ETM+ NDVI data derivatives are also expected to depict overall lower calculated
NDVI values as compared to Landsat 8 OLI equivalents (Roy et al., 2016). Although it is established
knowledge that atmospheric effects can reduce TOA product NDVI calculated values over land surface
areas (Roy et al., 2016), the alternative atmospherically corrected Surface Reflectance (SR) products may
not be appropriate for timeseries analyses over aquatic systems (Maciel et al., 2023). As the study area
included tidal wetland zones, TOA timeseries data were selected to minimize undesirable atmospheric
correction algorithm effects over coastal and inland water bodies. Consequently, modeled trend
analyses of the 30-year span of NDVI timeseries data needed take into consideration data variability
potentially attributable to sensor differences, atmospheric effects, and noisy data pixels.
Multi-year calculated NDVI averages for the Evergreen Forest class across non-overlapping sensor years
resulted in an estimated 0.05 overall NDVI value increase within Landsat 8 OLI scenes. This variation
seems in line with previously published studies (Yinghai Ke et al., 2015) and (Li P et al., 2014) that report
estimated deviations of ±0.05.
Image metadata file parameters were used to filter image scenes (Table App-1). The metadata filter
parameter 'CLOUD_COVER_LAND' was used to remove scene elements whose value exceeded 99%. A
secondary metadata filter parameter of 'DATA_TYPE' was used to include only 'L1TP' categorized data. A
51
-------
EPA 600/R-25/337 I April 2025
total of 47 image scenes qualified for removal from the timeseries due to the 'CLOUD_COVER_LAND'
parameter criteria and a total of 3 image scenes were excluded due the 'DATA_TYPE' parameter criteria.
Table App-1. Summarized totals for geographically corresponding unfiltered and filtered Landsat scenes for each Landsat
platform (with GEE Catalog identifier).
GEE Catalog Collection Identifier
Total # Unfiltered
Scenes
Total # Filtered
Scenes
Landsat 8 OLI
' LAN DSAT/LC08/C01/T1_RT_TOA'
516
500
Landsat 7 ETM+
'LANDSAT/LE07/C01/T1_RT_TOA'
1155
1134
Landsat 5 TM
'LANDSAT/LT05/C01/T1_TOA'
1256
1244
Landsat 4 TM
' LAN DSAT/LT04/C01/T1_TOA'
11
10
Data compilation into monthly NDVI products would generate an expected 372 total monthly composite
data layers for the study timeframe; however, gaps in data, attributable to satellite collection and/or
implemented data filters, reduced the total to 351 compiled monthly data layers. An initial filter was
used to identify image scenes that were flagged by their metadata to have been determined to have >=
99% cloud over land areas. These scenes, in their entirety, were removed from consideration in monthly
max NDVI calculations due to the inherent degree of reported cloud cover and the overall higher
likelihood of introducing undesirable data outliers to the timeseries.
The Landsat data quality assessment band ('BQA') for each image scene was used to generate and apply
a binary mask to remove cloud and cloud shadow pixels for each image scene in the timeseries
sequence. A cloud confidence bit value of 2, corresponding to a medium confidence level, was chosen to
generate image scene-specific pixel exclusion masks. Masked data pixels were assigned a flag value of -2
to facilitate monthly NDVI max value compositing and to facilitate post-identification of any remaining
flagged pixels. Data were sorted by timestamps and compiled to a monthly NDVI maximum value output
product. Proc Expand in SAS 9.2 (SAS®, Cary, North Carolina) with method "Join" was used to substitute
for missing values. The join function integrates linear interpolation lines of non-missing values over all
NDVI values. Interpolated values are then used to fill the missing. Data were subset to the study region
boundary and background (non-study) data areas were assigned a background flag value of -3. Data
were projected in GEE to European Petroleum Survey Group (EPSG): 3857 map projection and exported
in geoTiff raster file format at 30-meter pixel resolution for subsequent use.
References
Li, P., Jiang L., Feng, Z. 2014. Cross-Comparison of Vegetation Indices Derived from Landsat-7 Enhanced
Thematic Mapper Plus (ETM+) and Landsat-8 Operational Land Imager (OLI) Sensors. Remote
Sensing; 6(1): 310-329.
Maciel, D.A., Pahlevan, N., Barbosa, C.C.F., de Moraes de Novo, E.M.L., Paulino, R.S., Martins, V.S.,
Vermote, E., Crawford, C.J. 2023. Validity of the Landsat surface reflectance archive for aquatic
science: Implications for cloud-based analysis. Limnology and Oceanography Letters, 896: 850-858.
52
-------
PRESORTED STANDARD
POSTAGE & FEES PAID
EPA
PERMIT NO. G-35
Office of Research and Development (8101R)
Washington, DC 20460
Official Business
PenaltyforPrivate Use
$300
oEPA
United States
Environmental Protection
Agency
EPA 600/R-25/337 I April 2025 I
Recycled/Recyclable Printed on paper that contains a minimum of
50% postconsumer fiber content processed chlorine free
------- |