U.S. Department of the Interior

U.S. Geological Survey

Scientific Investigations Report 2015–5185

Prepared in cooperation with the North Dakota State Water Commission

Stochastic Model for Simulating Souris River Basin 

Precipitation, Evapotranspiration, and Natural 

Streamflow
Covers. Floodwaters of the Souris River surrounding homes (front) and railroad tracks (back) in 

Minot, North Dakota, June 27, 2011 (photograph by Brent R. Hanson, U.S. Geological Survey).
Stochastic Model for Simulating 

Souris River Basin Precipitation, 

Evapotranspiration, and Natural 

Streamflow

By Kelsey A. Kolars, Aldo V. Vecchia, and Karen R. Ryberg

Prepared in cooperation with the North Dakota State Water Commission

Scientific Investigations Report 2015–5185

U.S. Department of the Interior

U.S. Geological Survey
U.S. Department of the Interior

SALLY JEWELL, Secretary

U.S. Geological Survey

Suzette M. Kimball, Director

U.S. Geological Survey, Reston, Virginia: 2016

For more information on the USGS—the Federal source for science about the Earth, its natural and living 

resources, natural hazards, and the environment—visit http://www.usgs.gov or call 1–888–ASK–USGS.

For an overview of USGS information products, including maps, imagery, and publications, 

visit http://www.usgs.gov/pubprod/.

Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the 

U.S. Government.

Although this information product, for the most part, is in the public domain, it also may contain copyrighted materials 

as noted in the text. Permission to reproduce copyrighted items must be secured from the copyright owner.

Suggested citation:

Kolars, K.A., Vecchia, A.V., and Ryberg, K.R., 2016, Stochastic model for simulating Souris River Basin precipitation, 

evapotranspiration, and natural streamflow: U.S. Geological Survey Scientific Investigations Report 2015–5185, 55 p., 

http://dx.doi.org/10.3133/sir20155185.

ISSN 2328-0328 (online)
iii

Acknowledgments

We gratefully acknowledge Tim Fay and Mitch Weier of the North Dakota State Water Commission for their helpful reviews and Chanel Mueller of the U.S. Army Corps of Engineers, St. Paul 

District, for providing the natural (unregulated) streamflow data.

 We also wish to acknowledge U.S. Geological Survey employees Don Rosenberry and Stacey 

Archfield for their helpful reviews and Tara Gross for providing technical support.

v

Contents

Acknowledgments........................................................................................................................................iii

Abstract...........................................................................................................................................................1

Introduction.....................................................................................................................................................1

 Purpose and Scope .............................................................................................................................3

Analysis of Long-Term Climate Variability .................................................................................................4

Data and Methods ................................................................................................................................4

Climate Variability Results...................................................................................................................5

Stochastic Climate Model for Simulation of Precipitation, Temperature, and Potential 

Evapotranspiration...........................................................................................................................9

Data and Methods ................................................................................................................................9

Stochastic Simulation Results..........................................................................................................18

Water-Balance Model for Estimating Natural Streamflow...................................................................22

Input Data for Water-Balance Model..............................................................................................22

Streamflow..................................................................................................................................22

Watershed Delineation .............................................................................................................22

Soil Characteristics ...................................................................................................................24

Weather Data .............................................................................................................................24

 Water-Balance Model Description and Calibration.....................................................................26

Natural (Unregulated) Streamflow Routing...........................................................................30

Stochastic Natural Streamflow Model ....................................................................................................38

Methods................................................................................................................................................38

Stochastic Simulation Results..........................................................................................................39

Summary........................................................................................................................................................45

References Cited..........................................................................................................................................46

Appendix. Water-Balance Model Equations .......................................................................................51
vi

Figures

1. Map showing Souris River Basin showing locations of major reservoirs and

streamflow-gaging stations used for developing the stochastic natural streamflow

model ..............................................................................................................................................2

2. Graph showing annual peak streamflow for the Souris River above Minot, North

Dakota streamflow-gaging station 1904–2014 .........................................................................3

3. Map showing study area boundary, meteorological stations, and potential tree-ring

sites used for long-term climate analysis.................................................................................5

4. Map showing the five clusters of meteorological stations in the study area ....................6

5. Graphs showing multiresolution decomposition of tree-ring chronologies used in

the climate analysis......................................................................................................................7

6. Graph showing modeled and observed 12-year moving average precipitation for

cluster 1, season 2........................................................................................................................8

7. Graph showing modeled and observed 12-year moving average precipitation for

cluster 5, season 3........................................................................................................................9

8. Map showing locations of meteorological stations in and near the Souris River

Basin, including station groups used to develop the stochastic climate simulation

model for precipitation, temperature, and potential evapotranspiration and

additional stations used for calibrating the water-balance model ....................................10

9. Graphs showing time series plots of principal components for potential

evapotranspiration and precipitation for season 1, 1912–2011...........................................13

10. Graphs showing time series plots of principal components for potential

evapotranspiration and precipitation for season 2, 1912–2011...........................................14

11. Graphs showing time series plots of principal components for potential

evapotranspiration and precipitation for season 3, 1912–2011...........................................15

12. Graphs showing long-term simulated 12-year moving average precipitation for

season 2 and station group 1....................................................................................................21

13. Graphs showing long-term simulated 12-year moving average precipitation for

season 3 and station group 4....................................................................................................23

14. Map showing streamflow-gaging stations and corresponding subbasins,

designated by station identifier, used for developing the water-balance model.............25

15. Maps showing soil permeability and available water storage capacity ..........................28

16. Schematic of water-balance model for a generic grid cell.................................................29

17. Graph showing simulated and natural monthly streamflow for Long Creek at the

western crossing of the international boundary ...................................................................31

18. Graphs showing simulated and natural monthly streamflow for Long Creek at the

western crossing of the international boundary ...................................................................31

19. Graphs showing means of simulated and natural monthly streamflow............................33

20. Graphs showing standard deviations of simulated and natural monthly streamflow ....34

21. Graphs showing means of simulated and natural 10-day streamflow ..............................36

22. Graphs showing standard deviations of simulated and natural 10-day streamflow.......37

23. Graphs showing natural 10-day mean streamflow for 1946–2011 for the Souris River

below Rafferty Reservoir, Saskatchewan, Souris River above Minot, North Dakota,

and Wintering River near Karlsruhe, North Dakota, streamflow-gaging stations as

compared to the median, 90th percentile, and maximum recorded values of annual

maximum 10-day streamflow....................................................................................................40
vii

24. Graphs showing simulated 10-day mean streamflow for 100-year simulation period for

the Souris River below Rafferty Reservoir, Saskatchewan, Souris River above

Minot, North Dakota, and Wintering River near Karlsruhe, North Dakota,

streamflow-gaging stations as compared to the median, 90th percentile, and

maximum recorded values of annual maximum 10-day streamflow..................................41

25. Graphs showing simulated 10-day mean streamflow for 100-year simulation period

for the Souris River below Rafferty Reservoir, Saskatchewan, Souris River above

Minot, North Dakota, and Wintering River near Karlsruhe, North Dakota,

streamflow-gaging stations as compared to the median, 90th percentile, and

maximum recorded values of annual maximum 10-day streamflow..................................42

26. Graphs showing equilibrium frequency distributions for the Souris River below

Rafferty Reservoir, Saskatchewan, streamflow-gaging station for wet and dry

climate states ..............................................................................................................................43

27. Graphs showing equilibrium frequency distributions for the Souris River above

Minot, North Dakota, streamflow-gaging station for wet and dry climate states ...........44

28. Graphs showing equilibrium frequency distributions for the Souris River near

Westhope, North Dakota, streamflow-gaging station for wet and dry climate states ...45

Tables

1. Transformation parameters for time series variables used in the stochastic climate

simulation model.........................................................................................................................12

2. Results of principal components analysis for time series variables used in the

stochastic climate simulation model.......................................................................................12

3. Fitted multivariate periodic autoregressive models for adjusted principal

component time series variables. ..........................................................................................17

4. Statistics of simulated and observed precipitation and moisture deficit for dry 

period............................................................................................................................................19

5. Statistics of simulated and observed precipitation and moisture deficit for wet

period............................................................................................................................................20

6. Selected United States and Canadian streamflow-gaging stations used for

calibration of the water-balance model..................................................................................24

7. Selected United States and Canadian weather stations used for calibration of the

water-balance model .................................................................................................................26

8. Parameters from calibration of the monthly water-balance model using

reconstructed natural monthly streamflow from the U.S. Army Corps of Engineers......27

9. Lags applied to monthly streamflow routing model..............................................................32

10. Lags applied to 10-day streamflow routing model................................................................35
viii

Conversion Factors

Inch/Pound to International System of Units

Multiply By To obtain

Volume

acre-foot (acre-ft) 1,233 cubic meter (m3

)

acre-foot (acre-ft) 0.001233 cubic hectometer (hm3

) 

International System of Units to Inch/Pound

Multiply By To obtain

Length

centimeter (cm) 0.3937 inch (in.)

millimeter (mm) 0.03937 inch (in.)

kilometer (km) 0.6214 mile (mi)

kilometer (km) 0.5400 mile, nautical (nmi) 

Area

square kilometer (km2

) 247.1 acre

square kilometer (km2

) 0.3861 square mile (mi2)

Volume

cubic meter (m3

) 6.290 barrel (petroleum, 1 barrel = 42 

gal)

cubic meter (m3

) 264.2 gallon (gal) 

cubic meter (m3

) 0.0002642 million gallons (Mgal) 

cubic meter (m3

) 35.31 cubic foot (ft3)

cubic meter (m3

) 1.308 cubic yard (yd3) 

cubic meter (m3

) 0.0008107 acre-foot (acre-ft) 

Flow rate

cubic meter per second (m3

/s) 70.07 acre-foot per day (acre-ft/d) 

cubic meter per second (m3

/s) 35.31 cubic foot per second (ft3/s)

cubic meter per second (m3

/s) 22.83 million gallons per day (Mgal/d) 

centimeter per hour (cm/h) 0.3937 inch per hour (in/h)

Temperature in degrees Celsius (°C) may be converted to degrees Fahrenheit (°F) as °F = (1.8 × 

°C) + 32.

Datums

Horizontal coordinate information is referenced to the North American Datum of 1983 (NAD 83) 

and projected to the USA Contiguous Albers Equal Area Conic U.S. Geological Survey version.
Stochastic Model for Simulating Souris River Basin 

Precipitation, Evapotranspiration, and Natural Streamflow

By Kelsey A. Kolars, Aldo V. Vecchia, and Karen R. Ryberg

Abstract

The Souris River Basin is a 61,000-square-kilometer 

basin in the Provinces of Saskatchewan and Manitoba and 

the State of North Dakota. In May and June of 2011, recordsetting rains were seen in the headwater areas of the basin. 

Emergency spillways of major reservoirs were discharging at 

full or nearly full capacity, and extensive flooding was seen 

in numerous downstream communities. To determine the 

probability of future extreme floods and droughts, the U.S. 

Geological Survey, in cooperation with the North Dakota State 

Water Commission, developed a stochastic model for simulating Souris River Basin precipitation, evapotranspiration, and 

natural (unregulated) streamflow. Simulations from the model 

can be used in future studies to simulate regulated streamflow, 

design levees, and other structures; and to complete economic 

cost/benefit analyses. 

Long-term climatic variability was analyzed using treering chronologies to hindcast precipitation to the early 1700s 

and compare recent wet and dry conditions to earlier extreme 

conditions. The extended precipitation record was consistent 

with findings from the Devils Lake and Red River of the North 

Basins (southeast of the Souris River Basin), supporting the 

idea that regional climatic patterns for many centuries have 

consisted of alternating wet and dry climate states.

A stochastic climate simulation model for precipitation, 

temperature, and potential evapotranspiration for the Souris 

River Basin was developed using recorded meteorological 

data and extended precipitation records provided through treering analysis. A significant climate transition was seen around 

1970, with 1912–69 representing a dry climate state and 1970–

2011 representing a wet climate state. Although there were 

some distinct subpatterns within the basin, the predominant 

differences between the two states were higher spring through 

early fall precipitation and higher spring potential evapotranspiration for the wet compared to the dry state. 

A water-balance model was developed for simulating 

monthly natural (unregulated) mean streamflow based on 

precipitation, temperature, and potential evapotranspiration 

at select streamflow-gaging stations. The model was calibrated using streamflow data from the U.S. Geological Survey 

and Environment Canada, along with natural (unregulated) 

streamflow data from the U.S. Army Corps of Engineers. Correlation coefficients between simulated and natural (unregulated) flows generally were high (greater than 0.8), and the 

seasonal means and standard deviations of the simulated flows 

closely matched the means and standard deviations of the 

natural (unregulated) flows. After calibrating the model for 

a monthly time step, monthly streamflow for each subbasin 

was disaggregated into three values per month, or an approximately 10-day time step, and a separate routing model was 

developed for simulating 10-day streamflow for downstream 

gages.

The stochastic climate simulation model for precipitation, 

temperature, and potential evapotranspiration was combined 

with the water-balance model to simulate potential future 

sequences of 10-day mean streamflow for each of the streamflow-gaging station locations. Flood risk, as determined by 

equilibrium flow-frequency distributions for the dry (1912–69) 

and wet (1970–2011) climate states, was considerably higher 

for the wet state compared to the dry state. Future flood risk 

will remain high until the wet climate state ends, and for 

several years after that, because there may be a long lag-time 

between the return of drier conditions and the onset of a lower 

soil-moisture storage equilibrium. 

Introduction

The Souris River Basin is a 61,000-square-kilometer 

(km2) basin in the Provinces of Saskatchewan and Manitoba 

and the State of North Dakota (fig. 1). The basin topography 

consists of gently rolling prairie landscape and flat valleys. The main channel of the Souris River is approximately 

700 kilometers (km) in length, originating in southeast 

Saskatchewan, flowing southeast into North Dakota, and then 

continuing northward into Manitoba. The flow upstream from 

Minot, N. Dak., is affected by regulation from the Rafferty, 

Alameda, and Boundary Reservoirs in Saskatchewan and Lake 

Darling in North Dakota (fig. 1), which are operated in accordance with an international agreement between Canada and 

the United States (International Joint Commission, 1989). 

 Widespread record-setting flooding took place in the 

Souris River Basin from May to August of 2011. The flooding 
2 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

!

!

!!

!

!

!!

!

!

!

!

!

!

!

!

!!

!

#

#

#

#

Des Lacs

NWR Upper

Souris

NWR

J. Clark Salyer

NWR

#

!

(

(

((

(

(

((

(

(

(

(

(

(

(

(

((

(

*

*

*

*

Ç

Ç

Ç

Ç

Ç

Ç

Ç Ç

Ç

Ç

Ç

Ç Ç

Ç

Ç

Des Lacs

River

Wawanesa

Souris

Melita

Oxbow Estevan

Noonan Sherwood

Westhope

Willow City

Bantry

Minot

Foxholm

Karlsruhe

Verendrye

Upham

05NG001

05NG021 05NG007

05NF001 05ND004

05NB036

05NB001

05NA003 05124000

05123510 05123400

05122000

05120500

05120000 05117500

05116000

05116500

05114000 05113600

*

(

Deep River

Willow

River

Wintering

Assiniboine

Pipestone Creek

Moose

Souris River

Moose Jaw River

Long Creek

Mountain

Creek

Plum Creek

Lake Sakakawea

Creek

River

SASKATCHEWAN

MANITOBA

CANADA

UNITED STATES

EXPLANATION

MONTANA

NORTH DAKOTA

99° 100° 101° 102° 103° 104° 105°

50°49°48°

0 30 15 45 60 MILES0 50 25 75 100 KILOMETERS

National Wildlife RefugeReservoir and nameStreamflow-gaging station and identifier

05114000

Lake Winnipeg

AssiniboineRiver Basin

Red River of the North Basin

Souris  River  Basin

Devils Lake Basin

Alameda

Alameda

Lake Darling

Boundary

Rafferty

Figure 1. Souris River Basin showing locations of major reservoirs and streamflow-gaging stations used for developing the stochastic natural streamflow model.
Introduction 3

resulted from an unusually wet fall season in 2010, high 

snowfall during the subsequent winter season, and a series of 

record-setting rains during May and June 2011 in the headwater areas of the Souris River Basin in Saskatchewan and North 

Dakota (U.S. Department of Commerce, 2012). Emergency 

spillways of the upstream reservoirs were discharging at or 

near capacity, which resulted in extensive flooding to numerous downstream communities, including Minot, N. Dak., 

that had devastating damages with more than 4,000 homes 

flooded (U.S. Department of Commerce, 2012). Annual 

peak streamflow for the Souris River above Minot, N. Dak., 

(streamflow-gaging station 05117500) during 2011 was more 

than 736 cubic meters per second (m3/s), which was more than 

twice that of any previously recorded peak streamflow since 

1904 and more than five times the estimated 100-year postregulation peak streamflow (142 m3

/s; fig. 2; U.S. Geological 

Survey, 2014b). 

As a result of the record 2011 flood, a Souris River Basin 

Task Force convened to develop a plan of study for determining the likelihood of future severe flooding events, reviewing 

existing operating plans for upstream reservoirs, and improving downstream levees and other structures for mitigating 

future damages (International Souris River Board, 2013). One 

of the recommendations of the task force was to develop a stochastic model for simulating future streamflow and evaluating 

the probability of future severe floods or droughts. Streamflow 

modeling in the Souris River Basin is particularly challenging 

because of the extreme climatic variability and flat topography 

of the basin. As a result of the task force recommendation, 

the U.S. Geological Survey (USGS), in cooperation with the 

North Dakota State Water Commission, developed a stochastic 

model for simulating future Souris River Basin precipitation, 

evapotranspiration, and natural (unregulated) streamflow. 

Simulations from the model can be used in future studies to 

simulate regulated streamflow, design levees, and other structures, and complete economic cost/benefit analyses.

 Purpose and Scope

The purpose of this report is to describe a stochastic natural streamflow model and present results of model simulations 

for evaluating the probability of future floods in the Souris 

River Basin. The stochastic natural streamflow model was 

developed through a series of four stages that are described in 

subsequent sections of this report:

1. analysis of long-term climatic variability in the Souris 

River Basin and surrounding region,

2. development of a stochastic climate model for simulating 

climatic inputs,

3. development of a water-balance model for simulating 

natural (unregulated) monthly streamflow in response to 

climatic inputs, and 

4. development of a natural (unregulated) streamflow routing model for simulating streamflow at a three-per-month 

(approximately 10-day) time step.

0

100

200

300

400

500

600

700

800

1900 1910 1920 1930 1940 1950 1960

Year

1970 1980 1990 2000 2010 2020

Annual peak streamflow, in cubic meters per second

Figure 2. Annual peak streamflow for the Souris River above Minot, North Dakota streamflow-gaging station (05117500) 

1904–2014.
4 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

Analysis of Long-Term Climate 

Variability

Long-term climate variability was analyzed using treering data to extend historical precipitation records from the 

Souris River Basin and surrounding areas back to the early 

1700s and to compare recent (since the late 1800s) “extreme” 

wet and dry conditions to earlier conditions. The methods and 

results summarized in this section are described in detail in 

Ryberg (2015). A brief overview is given in this report to provide background information for the subsequent sections. 

Data and Methods

To evaluate climate variability for the Souris River Basin 

in a larger regional context, historical precipitation and temperature data were obtained from United States and Canadian 

meteorological stations in an area extending beyond the Souris 

River Basin boundary and including parts of other basins in the 

region, including the Missouri River, Red River of the North, 

Assiniboine River, and Saskatchewan River (fig. 3). Monthly 

precipitation and temperature (average, maximum, and minimum) data were downloaded for Canadian stations from the 

Adjusted and Homogenized Canadian Climate Data (AHCCD; 

Environment Canada, 2014) and were downloaded for U.S. stations from the United States Historical Climatology Network 

(USHCN; Menne and others, 2014). The AHCCD and USHCN 

datasets are designed to minimize potential trends because of 

changes in instrumentation, station locations, and other factors 

unrelated to natural climate variability or climate change.

Tree-ring data were downloaded from The International 

Tree-Ring Data Bank (ITRDB; National Oceanic and Atmospheric Administration National Climate Center, 2014). As 

described in Ryberg (2015), after careful analysis of data 

from sites within and near the study boundary, four sites were 

selected for inclusion in the analysis (fig. 3): Boundary Bog in 

Saskatchewan, Theodore Roosevelt National Park and Burning 

Coal Vein in North Dakota, and Cedar Butte in South Dakota.

Monthly precipitation and temperature data for each of the 

AHCCD and USHCN meteorological stations were aggregated 

to obtain seasonal total precipitation and average temperature 

for three 4-month seasons: season 1 (November–February), 

season 2 (March–June), and season 3 (July–October). These 

seasons were selected because climatic conditions in each 

season have a distinct and important effect on runoff for the 

Souris River Basin. Precipitation during season 1 generally 

is trapped in frozen soils or stored in the basin as snowpack 

and either replenishes soil moisture deficits from the previous 

summer or provides snowmelt runoff. Season 2 generally is the 

period when peak streamflow from snowmelt runoff and (or) 

spring rainfall takes place, and season 3 generally is a period of 

low flows and moisture deficits as potential evapotranspiration 

(PET) exceeds precipitation and soil moisture is depleted. 

Hierarchical agglomerative cluster analysis (HACA; 

Maechler and others, 2014) was used to identify five distinct 

climatic clusters of sites based on mean seasonal precipitation (1900–2010; fig. 4). Precipitation and temperature were 

initially considered for the cluster analysis, but precipitation 

was by far the dominant of the two variables; thus, temperature was not used (Ryberg, 2015). The position of the Souris 

River Basin and its partitioning into multiple clusters (fig. 4) 

highlights the complexity and variability of climatic conditions in the basin despite its small size. 

To focus on regional-scale temporal variability, precipitation data for each season and each cluster were averaged 

over the stations in the cluster to obtain 15 (5 clusters times 3 

seasons) yearly time series for the period 1888 through 2010; 

furthermore, to remove high-frequency variability and focus 

on long-term (multidecadal) variability, the yearly time series 

were averaged over 12-year moving windows. A 12-year 

window was used because a number of studies have identified 

quasi-periodic signals with wavelengths of 12 years or less 

for precipitation in the north-central United States and southcentral Canada during various seasons (Ault and St. George, 

2010; Small and Islam, 2008, 2009; Garbrecht and Rossel, 

2002; Yang and others, 2007). 

To hindcast the seasonal precipitation, each of the 15 

smoothed yearly time series described in the previous paragraph were modeled based on the annual tree-ring chronologies. Multiresolution decomposition (Bunn and others, 2014; 

Whitcher, 2013) was used to express the annual tree-ring data 

for each site in terms of additive frequency components or 

wavelet voices (fig. 5). Because the precipitation data were 

smoothed using a 12-year window, the high-frequency wavelet 

voices (D1, D2, and D3, fig. 5), corresponding to wavelengths 

of 2, 4, and 8 years, were not used in the analysis. Also, 

because the tree-ring chronologies were not long enough to 

have confidence in the D7 wavelength (128 years), D7 was not 

used in the analysis; therefore, for all four tree-ring chronologies, three wavelet voices (D4, D5, and D6), corresponding 

to wavelengths of 16, 32, and 64 years, were used to estimate 

precipitation. Thus, there were 12 potential explanatory variables (plus an intercept) used to model each of the 15 yearly 

precipitation time series. The calibration period (the time 

interval for which both tree-ring and precipitation data were 

available) extended from 1888, 1889, or 1890 (depending on 

the precipitation group) to 1990, the year shortly after which 

most of the tree-ring chronologies ended. The beginning of the 

hindcast period (the period during which only tree-ring data 

were available) was chosen to be 1700. Although three of the 

tree-ring chronologies extended well back into the 1600s, the 

Boundary Bog chronology, which began about 1700, turned 

out to be the most frequently used explanatory variable. 

As described in Ryberg (2015), all-subsets regression 

along with extensive model selection and verification criteria 

were used to select the best models for explaining variability 

in each of the precipitation time series based on the tree-ring 

wavelet voices. To reduce serial correlation, the regression 
Analysis of Long-Term Climate Variability 5

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

)

)

)

)

Cedar Butte

Boundary Bog

Burning Coal Vein

Theodore 

Roosevelt 

National Park

Qu'Appelle River

Missouri River

River

Souris

MANITOBA

SASKATCHEWAN

MONTANA

WYOMING

MINNESOTA

SOUTH DAKOTA

NORTH DAKOTA Red River of the North

Saskatchewan River

Assiniboine

River

120° 110° 100° 90°

50°

40°

EXPLANATION

) Tree-ring sites

! Meteorological stations

Souris River Basin

Four-degree buffer

Base map from Commission for Environmental Cooperation North American Environmental Atlas digital data, 2010

North American Datum of 1983

Universal Transverse Mercator projection, zone 14

CANADA

UNITED STATES

Figure 3. Study area boundary (four-degree buffer around Souris River Basin), meteorological stations, and 

potential tree-ring sites used for long-term climate analysis. Modified from Ryberg (2015).

models were based on thinned (every fourth year) values of 

the annual time series.

Climate Variability Results

The regression analysis for estimating the smoothed 

yearly precipitation data for each season and cluster based 

on the tree-ring data indicated that, with few exceptions, the 

tree-ring wavelet voices for the selected sites explained most 

of the multidecadal variability in seasonal precipitation (table 

5.1 in Ryberg, 2015). The regression models explained more 

than 60 percent of the variability in precipitation (based on 

the adjusted coefficient of determination) for cluster 1 (all 

3 seasons), cluster 2 (seasons 1 and 3), cluster 3 (season 1), 

cluster 4 (seasons 1 and 3), and cluster 5 (all 3 seasons); and 

many of these models explained more than 80 percent of the 

variability. Of the remaining models, the models for cluster 2 

(season 2) and cluster 4 (season 2) explained 58 and 54 percent, respectively, of the variability. Only cluster 3 (seasons 

2 and 3) had a relatively low coefficient of determination (32 

and 45 percent, respectively), and these were the only models 

that failed the regression diagnostic tests; therefore, with the 

exception of cluster 3 (seasons 2 and 3), the tree-ring chronologies were good predictors of seasonal precipitation for all 

five clusters. The most frequently used explanatory variables 

included wavelet voice D6 for Burning Coal Vein (10 models), 
6 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

110° 105° 100° 95°

55°

50°

45°

" Meteorological stations in 

 cluster analysis

Climate-based clusters Souris River Basin

1 Southeast

2 South-central

3 Southwest

4 Northwest

5 Northeast

Base map from Commission for Environmental Cooperation North American Environmental Atlas digital data, 2010

North American Datum of 1983

Universal Transverse Mercator projection, zone 14

EXPLANATION

MANITOBA SASKATCHEWAN

MONTANA

WYOMING

MINNESOTA

SOUTH DAKOTA

NORTH DAKOTA

CANADA

UNITED STATES

Figure 4. The five clusters of meteorological stations in the study area. Modified from Ryberg (2015).
Analysis of Long-Term Climate Variability 7

Multiresolution decomposition of Boundary Bog

D7

D6

D5

D4

D3

D2

D1

Multiresolution decomposition of Theodore Roosevelt National Park

Multiresolution decomposition of Burning Coal Vein

Multiresolution decomposition of Cedar Butte

1600 1700 1800

Year

1900 2000

D7 128

D6 64

D5 32

D4 16

D3 8

D2 4

D1

D7

D6

D5

D4

D3

D2

D1

D7

D6

D5

D4

D3

D2

D1

2

128

64

32

16

8

128

64

32

16

8

128

64

32

16

8

4

2

4

4

2

2

Years

Wavelet detail vectors

Figure 5. Multiresolution decomposition of tree-ring chronologies used in the climate analysis. Modified from Ryberg 

(2015).
8 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

Boundary Bog (8 models), and Theodore Roosevelt National 

Park (8 models); and wavelet voice D5 for Cedar Butte (9 

models), Boundary Bog (8 models), and Theodore Roosevelt 

National Park (6 models). 

The extended precipitation time series for cluster 1 

(southeast) and season 2 (March–June) is shown in figure 6 

along with the observed precipitation data. The coefficient 

of determination for this model was 62 percent (based on the 

thinned data). Observed season 2 precipitation for this group 

was highest at the end of the record (about 1993 to present) 

and during about 1890–1910. Similar (in magnitude) but 

longer duration wet periods were seen about 1820–50 and 

1740–70. There were two extreme (compared to the observed 

record) season 2 droughts during the 1860s and early 1700s. 

These extreme wet and dry periods prior to the precipitation 

record were confirmed in Ryberg (2015) based on historical 

information and comparison with other paleoclimatic studies 

from the region.

The extended precipitation time series for cluster 5 

(northeast) and season 3 (July–October) is shown in figure 7 

along with the observed 12-year moving average precipitation 

data. The coefficient of determination for this model was 

81 percent. Observed 12-year moving average precipitation 

for this season is extremely high from about 1980 to present, remaining well above previously recorded precipitation. 

A shift toward higher summer and fall precipitation in about 

1980 has been documented in a number of previous studies 

in the southeastern part of the study area and elsewhere in the 

north-central United States and south-central Canada (Garbrecht and Rossel, 2002; Shapely and others, 2005; Vecchia, 

2008; Small and Islam, 2008). Based on figure 7, this recent 

wet summer/fall period extends well into southern Manitoba 

and southeastern Saskatchewan, including a large part of the 

Souris River Basin (fig. 4). The extended precipitation record 

indicates a similar (though perhaps shorter duration) wet 

period during about 1830–50 and another extreme wet period 

during the early 1700s (fig. 7). Extreme dry summer/fall 

conditions are indicated during the 1740s and 1790s, and these 

droughts are similar to the drought of the 1930s (fig. 7). These 

extreme wet and dry periods, indicated by the extended precipitation data, were verified in Ryberg (2015) based on historical 

information and comparison with other paleoclimatic studies.

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●●●

●

●

●

●

●

●

●●

●

●

●

●

●

●●●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●

●

●

●

●●●

●

●

●

●

●

●●●

●

●

●

●

●

●

●●●●●

●

●

●

●

●

●

●

●

●

●

●●●●●●●●●●

●

●

●

●

●

●

●●●

●

●

●

●

●

●

●

●

●●●●●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●●●●●●●●●●●●●

●

●

●

●

●

●

●●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●●

●

●

●

●

●

●

●

●

●

●●●●

●

●●●●

●

●

●

●

●

●

●

●

●

●

●

●●

●

●

●

●

●

●

●

●●●

●

●

●

●

●

●●

●

●

●

●

●

●

●

●

●●●

●

●

●

●

●

●

●

●

●

●●●●

●

●

●

●

●

●

●

●

●●●

●

●

●

●

●

●

●

●

●

●●●

●

●

●

●

●

●●●●●●

●

●

●

●

●

●

●

●

●

●●

●

●

●

●

●

●

●

●

●

●●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●

●

●

●

●

●

●

●

●●

●

●

●

●

●

●

●

●

●

●

●

●

●●

●

●

●●

●

●●

●

●

●●

●●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●

●

●

●●

●

●

●●

●

●

●

●

●

●●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●●

●

180

160

200

220

240

260

Total seasonal precipitation, in millimeters

Midpoint of 12-year moving average

1675 1725 1775 1825 1875 1925 1975 2025

Observed precipitation Modeled precipitation

Figure 6. Modeled and observed 12-year moving average precipitation for cluster 1 (southeast, representing much 

of the Red River of the North Basin in the United States), season 2 (March-June).
Stochastic Climate Model for Simulation of Precipitation, Temperature, and Potential Evapotranspiration 9

180

200

220

240

260

Midpoint of 12-year moving average

Total seasonal precipitation, in millimeters

●●●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●●●●●●

●

●

●

●

●

●

●

●

●●

●

●

●

●

●

●

●

●

●

●

●

●●●

●

●

●

●

●

●

●

●●●●●●●●●●●●

●

●

●

●

●

●

●

●

●

●●●

●

●

●

●

●

●

●

●●●●●●●●●●●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●●●●●●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●●●●●

●

●

●

●

●

●

●

●

●

●

●

●●●

●

●

●

●

●

●

●

●

●●

●

●

●

●

●

●

●

●

●●●

●

●

●

●●●●●●●●●●●

●

●

●

●

●●●●

●

●

●

●

●

●

●

●

●●

●

●

●

●

●

●

●

●

●

●

●

●●●

●

●

●

●

●

●●●●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●●●●●●

●

●

●

1675 1725 1775 1825 1875 1925 1975 2025

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●

●●

●

●

●

●

●

●

●

●

●

●

●

●●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●

●

●

●

●

●

●

●

●

●

●

●●

●

●

●

●●

●

●

●

●●

●

●

●

●

●

●

●

●

●●

●

●

●

●

●

●

●

●

●●

●

●

●

●

Observed precipitation Modeled precipitation

Figure 7. Modeled and observed 12-year moving average precipitation for cluster 5 (northeast, representing 

southern Manitoba and southeastern Saskatchewan), season 3 (July–October).

Stochastic Climate Model for 

Simulation of Precipitation, 

Temperature, and Potential 

Evapotranspiration

The water-balance model for estimating natural streamflow, described in subsequent sections of this report, uses 

monthly precipitation, temperature, and PET to estimate runoff 

from selected headwater basins and intervening subbasins of 

the Souris River Basin. The water-balance model was calibrated and verified using historical climate inputs and streamflow data; however, the objective of this study is to generate 

potential future streamflow. Therefore, a stochastic climate 

model was developed to simulate future precipitation, temperature, and PET over the Souris River Basin that reproduces 

both the short-term variability in the historical meteorological 

records and the long-term variability in the extended tree-ring 

precipitation records. The stochastic climate simulation model 

is described in this section. 

Data and Methods

The meteorological data described previously in the section “Analysis of Long-term Climate Variability” were examined to select stations in and near the Souris River Basin with 

complete or nearly complete records of monthly precipitation 

and average temperature for 1912–2011. The starting year 

(1912) was selected because many of the stations had missing 

or partially missing records prior to that date. The large areal 

extent of the regions used previously (fig. 4) contained many 

stations and thus spatial averages over the regions could be 

estimated for earlier years despite missing records for some 

stations; however, for this analysis it was important to use 

stations in and closely surrounding the Souris River Basin, and 

the stations needed to have a common period of record. There 

were 27 stations closest to the basin that were initially considered (fig. 8). Of these stations, 16 had continuous records 

for 1912–2011. The 16 continuous-record sites were split into 

4 station groups of 4 stations each (fig. 8). Note that the station groups defined here differ from the clusters defined in the 

previous section of this report. This grouping of stations was 
10 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

k

k

k

k

!

!

!

!

!

!

!

!

!k!

MANITOBA

SASKATCHEWAN

MONTANA

NORTH DAKOTA

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

j

k

j

j

j

k

k

.

>

.

>

>

.

>

.

k

Mandan experiment station

100° 105°

51°50°49°48°47°

0 60 80 120 KILOMETERS 30

.j>k

" 0 40 60 MILES 20

Vida

Medicine Lake3 southeast

Crosby

Estevan

Regina Indian Head

Kelliher

Tonkin

Bangor

Manor

Pierson

Bottineau

Willow City

Ninette

Cypress River

Brandon

PortagePrairie

Langruth West

NeepawaMurray

Morden

Langdonexperimentalfarm

Towner2 northeast

Moosomin

Birtle

Yellow Grass

Ceylon

Northeast Group and station nameNorthwest Group and station nameSoutheast Group and station nameSouthwest Group and station nameWeather station

1-Degree buffer

Subbasin

EXPLANATION

CANADA

UNITED STATES

Figure 8. Locations of meteorological stations in and near the Souris River Basin, including station groups used to develop the stochastic climate simulation model for precipitation, temperature, and potential evapotranspiration (northeast, northwest, southeast, and southwest groups) and additional stations used for calibrating the water-balance model.
Stochastic Climate Model for Simulation of Precipitation, Temperature, and Potential Evapotranspiration 11

the most compact set of stations with continuous record that 

was sufficient to generate the climatic inputs required for the 

water-balance model. Some of these stations had occasional 

missing values but none had large (greater than 1 year) periods 

of missing data.

Monthly PET values were estimated using the Hamon 

method (Lu and others, 2005; Huffman and others, 2011) and 

modified as described in the following section “Water-Balance 

Model for Estimating Natural Streamflow” of this report. 

The Hamon method estimates PET as a function of monthly 

average temperature, position (latitude), and time of year 

(radiation as it relates to declination of the Earth). For a given 

spatial position and month, PET computed using this method 

is a one-to-one, monotonically increasing function of monthly 

average temperature; therefore, because monthly temperature 

can be computed directly from monthly PET, the stochastic 

climate simulation model was developed for precipitation and 

PET. Temperature, which is required for the snowmelt algorithm used in the water-balance model, was computed directly 

from PET as described later in the section “Stochastic Natural 

Streamflow Model” of this report.

The stochastic climate model for simulating monthly 

precipitation and PET for the 16 sites contained in the groups 

shown in figure 8 was developed in two stages. First, a model 

was developed for group-averaged seasonal data for each of 

the four station groups and the same three seasons used in the 

section “Analysis of Long-Term Climate Variability.” Then, 

the group-averaged seasonal data were disaggregated, as 

described later in the section “Stochastic Natural Streamflow 

Model,” to obtain monthly values for individual sites. The 

group-averaged seasonal data is designated as

PET y P y GS GS , , ( ); ( ) (1)

where

PET y G S, ( ) is the potential evapotranspiration for station 

group G, season S, and year y; P y G S, ( ) is the precipitation for station group G, season 

S, and year y;

G equals 1 (southeast), 2 (southwest), 3 

(northwest), and 4 (northeast);

S equals 1 (November–February), 2 (March–

June), and 3 (July–October); and

y equals 1912, 1913, ..., 2011.

Henceforth, the variables in equation 1 will be referred 

to as time series variables. The first step in the analysis was 

to transform each of these 24 time series variables (2 climate 

variables times 4 groups times 3 seasons) so that the transformed values had an approximate standard normal distribution. This was accomplished using power transformations:

ZPET

PET M

SD

ZP

P c M

SD

r r

= [ ] − − = c [ ] − − ; , (2)

where 

ZPET is the transformed value of PET;

ZP is the transformed value of P;

c and r are the transformation parameters;

M is the sample mean, values provided in table 

1; and

SD is the sample standard deviation, values 

provided in table 1. 

The transformation parameters (table 1) were selected 

as follows. For c=0, values r=1, 1/2, 1/3, 1/4, and 1/5 were 

applied, and the value for which the skewness coefficient of 

the transformed data was closest to zero (in absolute value) 

was selected. The resulting values were r=1 (no transformation) for PET for season 3 and r=1/3 (season 1) or r=1/2 (seasons 2 and 3) for P. For PET during seasons 1 and 2, the skewness coefficient was positive and did not approach zero as r

was decreased. In those cases, c=10 was used for season 1 and 

c=100 for season 2, and the value r=1/4 resulted in a skewness coefficient close to zero. The M and SD are the sample 

means and standard deviations of the transformed variables for 

1912–2011. Normal probability plots of ZPET and ZP were 

examined and in all cases the standard normal approximation 

was deemed adequate.

The next step in the model development was to use 

principal components analysis to evaluate spatial patterns in 

the transformed time series variables that could be used to 

simplify the model. For each season, principal components 

are linear combinations of the time series variables over the 

station groups; for example, for PET, the principal components 

can be expressed as

ZPET y a ZPET y j j S jG G S

G

,

*

, ( ) = = ( ), , , , =

∑1

4

 for 1 234 (3)

Corr ZPET y ZPET y j k j S k S { ( ), ( )} , ,

*

,

* = ≠ 0 (4)

Var ZPET y Var ZPET y

Var ZPET y Var ZPE

S S

S

{ ( )} { ( )}

{ ( )} {

,

*

,

*

,

*

1 2

3

>

> > T y 4,S

* ( )} (5)

Var ZPET y Var ZPET y j S

j

S

G

{ ( )} { ( )} ,

* = G,

= =

∑ ∑ 1

4

1

4

(6)

where

ZPET y j S,

* ( ) is the transformed value of PET for principal 

component j, season S, and year y;

ZPET y G S, ( ) is the transformed value of PET for station 

group G, season S, and year y; and

ajG is the coefficient for principal component j

and station group G.

The coefficients of the principal components (a in eq. 3) 

are computed such that the components are mutually uncorrelated (eq. 4), the variance of the components decreases 

as the component number (j) increases (eq. 5), and the total 

variance of the principal components equals the total variance 
12 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

Table 1. Transformation parameters for time series variables used in the stochastic climate simulation model.

[c and r, transformation parameters defined in equation 2; M, sample mean; SD, sample standard deviation; SE, southeast; SW, southwest; NW, northwest; NE, 

northeast]

Season Group Potential evapotranspiration Precipitation

c r M SD c r M SD

1

SE 10 0.25 2.24 0.13 0 0.33 3.61 0.47

SW 10 0.25 2.30 0.14 0 0.33 3.64 0.43

NW 10 0.25 2.10 0.16 0 0.33 4.50 0.41

NE 10 0.25 2.17 0.14 0 0.33 4.32 0.44

2

SE 100 0.25 3.26 0.12 0 0.5 13.48 2.01

SW 100 0.25 3.31 0.13 0 0.5 12.93 2.08

NW 100 0.25 3.14 0.13 0 0.5 13.63 2.17

NE 100 0.25 3.25 0.12 0 0.5 14.12 1.99

3

SE 0 1.0 301.7 17.5 0 0.5 13.66 1.90

SW 0 1.0 309.0 18.7 0 0.5 12.32 1.75

NW 0 1.0 274.5 16.9 0 0.5 13.60 1.89

NE 0 1.0 298.5 17.3 0 0.5 14.53 1.90

Table 2. Results of principal components analysis for time series variables used in the stochastic climate simulation model.

[PC, principal component; SE, southeast; SW, southwest; NW, northwest; NE, northeast]

Season PC 

number

Potential evapotranspiration Precipitation

Variance Coefficients for station groups Variance Coefficients for station groups

SE SW NW NE SE SW NW NE

1

1 3.84 0.50 0.50 0.50 0.50 2.71 0.48 0.51 0.50 0.51

2 0.07 0 -0.77 0.16 0.62 0.56 0.72 -0.22 -0.64 0.16

3 0.04 -0.64 0 0.75 -0.16 0.44 -0.19 -0.69 0.19 0.67

4 0.01 -0.58 0.40 -0.40 0.59 0.25 -0.46 0.47 -0.56 0.51

2

1 3.77 0.50 0.50 0.50 0.50 2.93 0.50 0.50 0.47 0.53

2 0.10 0.48 -0.47 -0.53 0.51 0.59 0.59 -0.34 -0.65 0.34

3 0.06 -0.30 -0.65 0.62 0.33 0.32 -0.11 -0.77 0.50 0.39

4 0.02 0.65 -0.34 0.30 -0.61 0.12 0.63 -0.22 0.32 -0.67

3

1 3.74 0.51 0.49 0.47 0.53 2.76 0.51 0.51 0.45 0.53

2 0.16 0.50 -0.41 -0.62 0.45 0.65 0.51 -0.27 -0.72 0.38

3 0.08 0 -0.75 0.61 0.22 0.34 0.23 -0.81 0.53 0.10

4 0.04 0.70 -0.13 0.15 -0.68 0.21 0.65 0.11 0 -0.75

of the original data (eq. 6). The statistical software R function “princomp” (R Core Team, 2015) was used to evaluate 

principal components for each season, and the results are 

given in table 2. For both ZPET and ZP and all three seasons, 

the first principal component (PC1) consisted of the average of all four station groups (multiplied by a constant). PC1 

explained almost all (about 95 percent) of the variability for 

ZPET and most the variability (from 67 to 73 percent) for 

ZP. The second principal component (PC2) was similar for 

all three seasons and both ZPET and ZP and consisted of a 

weighted average of the eastern groups minus a weighted 

average of the western groups, or an east-west contrast. PC2 

explained less than 4 percent of the variability for ZPET and 

from about 14 to 16 percent of the variability for ZP. The 

third principal component (PC3) also was similar for all three 

seasons and both ZPET and ZP and consisted of a weighted 
Stochastic Climate Model for Simulation of Precipitation, Temperature, and Potential Evapotranspiration 13

average of the northern groups minus a weighted average of 

the southern groups, or a north-south contrast. PC3 explained 

less than 2 percent of the variability for ZPET and from about 

8 to 11 percent of the variability for ZP. The fourth principal 

component (PC4) explained less than 1 percent of the variability for ZPET and less than 6 percent of the variability for ZP.

Based on the results presented in the previous paragraph, 

a model for the first three principal components for both ZPET

and ZP was developed. Although, for ZPET, the second and 

third components explain a small part of the variability, these 

components were included because they were similar (in terms 

of the coefficients) to those for ZP and were important for 

maintaining cross-correlation between ZPET and ZP. PC4 is 

not included in any of the subsequent plots. PC4 was assumed 

to consist of white noise (independently distributed normal 

random variables) with mean zero and variance given by the 

corresponding values from table 2. In the simulation model 

described later in this section, all four principal components 

needed to be generated in order to be able to back-compute 

ZPET and ZP from the principal components.

Time series plots of the first three principal components 

for each of the seasons are shown in figures 9–11. A stochastic model needed to be developed to jointly simulate these 

18 time series variables (in addition to the 6 variables for PC4, 

which were simulated as described in the previous paragraph). 

The model needed to reproduce the statistical properties of 

the data, including any important correlations between values 

from one season to the next or between values for PET and 

Significant step trend

1970 break point in step trend

EXPLANATION

Principal component 1 (average)

−6

−5

−4

−3

−2

−1

0

1

2

3

4

5

6

−4

−3

−2

−1

0

1

2

3

4

Principal component 2 (east-west contrast)

−1

0

1

−3

−2

−1

0

1

2

3

Principal component 3 (north-south contrast)

−1

0

1

1900 1920 1940 1960

Year Year

1980 2000 2020

−3

−2

−1

0

1

2

3

1900 1920 1940 1960 1980 2000 2020

Principal component value (dimensionless)

Potential evapotranspiration Precipitation

Figure 9. Time series plots of principal components for potential evapotranspiration and precipitation for season 1

(November–February), 1912–2011.
14 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

Principal component 1 (average)

Principal component 2 (east-west contrast)

−1

0

1

−3

−2

−1

0

1

2

3

Principal component 3 (north-south contrast)

−1

0

1

1900 1920 1940 1960 1980 2000 2020 1900 1920 1940 1960 1980 2000 2020

Principal component value (dimensionless)

Potential evapotranspiration Precipitation

−7

−6

−5

−4

−3

−2

−1

0

1

2

3

4

5

6

7

−5

−4

−3

−2

−1

0

1

2

3

4

5

−2

−1

0

1

2

Significant step trend

1970 break point in step trend

EXPLANATION

Year Year

Figure 10. Time series plots of principal components for potential evapotranspiration and precipitation for season 2 (March–

June), 1912–2011.
Stochastic Climate Model for Simulation of Precipitation, Temperature, and Potential Evapotranspiration 15

Principal component 1 (average)

Principal component 2 (east-west contrast)

Principal component 3 (north-south contrast)

1900 1920 1940 1960 1980 2000 2020 1900 1920 1940 1960 1980 2000 2020

Principal component value (dimensionless)

Potential evapotranspiration Precipitation

−5

−4

−3

−2

−1

0

1

2

3

4

5

−5

−4

−3

−2

−1

0

1

2

3

4

5

−2

−1

0

1

2

−2

−1

0

1

2

−1

0

1

−2

−1

0

1

2

Significant step trend

1970 break point in step trend

EXPLANATION

Year Year

Figure 11. Time series plots of principal components for potential evapotranspiration and precipitation for season 3 (July–

October), 1912–2011.
16 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

precipitation. The model also needed to reproduce any nonstationarity (such as trends in the mean) or long-range dependence (small deviations of the central tendency from zero 

that persist for long periods of time before reversing direction). Based on the climate analysis described in the section 

“Analysis of Long-Term Climate Variability” in this report, 

and on previous studies from the Devils Lake Basin in North 

Dakota (Bluemle, 1996; Wiche and others, 2000; Vecchia, 

2002, 2008; Hoerling and others, 2010), the Red River of the 

North Basin (Rannie, 1998; St George and Nielson, 2003), and 

elsewhere (St George and Neilson, 2002; Shapely and others, 

2005; Ryberg and others, 2014; Yonetani and Gordon, 2001), 

there seems to be strong support for long-range dependence; 

however, evidence for trends or permanent changes has not 

been firmly established for central North America (Milly and 

others, 2005; Hirsch and Ryberg, 2012; Ryberg and others, 

2014). Of particular relevance to this study is evidence from 

the Devils Lake Basin, adjacent to and east of the Souris 

River Basin (fig. 1), in which an abrupt and highly significant 

increase in precipitation during summer and fall was seen in 

the late 1970s along with an increase in spring and summer 

lake evaporation, which is highly correlated with PET (Vecchia, 2008). Based on a comparison between simulated lake 

levels and 5,000 years of reconstructed lake levels from paleo 

lake studies (Murphy and others, 1997), Vecchia (2008) concluded that climate in the Devils Lake Basin had undergone 

abrupt transitions from “dry” to “wet” (such as was seen in the 

late 1970s) and “wet” to “dry,” on average, about once every 

150 years. The duration and timing of the transition periods 

are highly variable and seemingly random.

Because of the results for the Devils Lake Basin in Vecchia (2008), it was postulated that the Souris River Basin may 

have undergone a similar abrupt climate transition at about 

the same time; therefore, for each of the time series variables 

in figures 9–11, two-sample t-tests (Helsel and Hirsch, 2002) 

were used to test for step trends in the mean at three transition 

times: 1960, 1970, and 1980. There were several significant 

step trends in each case, but the trends were most significant 

for the 1970 transition time. This differed somewhat from 

Devils Lake, in which 1980 produced the most significant step 

trends; however, in either case the step trends may have taken 

place sometime between 1970 and 1980 with the exact year 

subject to speculation. The significant (p-value, or obtained 

significance level, less than 0.01) step trends for the 1970 

transition time for this study are indicated by the red lines in 

figures 9–11. For PC1 (regional average of the four groups), 

there were significant uptrends in PET for season 1 (fig. 9), 

PET and precipitation for season 2 (fig. 10), and precipitation for season 3 (fig. 11). These trends are consistent with 

uptrends in precipitation and lake evaporation observed for the 

Devils Lake Basin and may be indicative of a large regional 

climate “shift” that was seen sometime during the early to late 

1970s. There also was a significant downward step trend in 

PC2 (east-west contrast) and a significant upward trend in PC3 

(north-south contrast) for PET for season 2 (fig. 10), indicating 

that the increase in PET was largest for the northern and western part of the region (that is, the NW group). For season 3 

PET, there was a significant downtrend for PC2, despite no 

trend in PC1, which indicates that, although average PET for 

all four groups did not change, there were small increases in 

the southwest and northwest groups offset by small decreases 

in the southeast and northeast groups. The implications of 

these changes in terms of original (untransformed) PET and 

precipitation and their impact on the water balance will be 

discussed later in this section. 

To simulate the principal component time series, the 

step changes described in the previous paragraph were subtracted out to obtain adjusted time series for each principal 

component: 





EyZPET y y

P y ZP y y

j S j S j S

j S j S j S

, ,

*

,

, ,

*

,

( ) ( ) ( );

( ) ( ) ( )

= −

= −

α

β

(7)

 j S D j S y W j S y y A I A I , ; , ; , ( ) = + 1912≤ ≤1969 1970≤ ≤2011 (8)

 j S D j S y W j S y y B I B I , ; , ; , ( ) = + 1912≤ ≤1969 1970≤ ≤2011 (9)

where 

E y 

j S, ( ) is the adjusted value of ZPET* for principal 

component j, season S, and year y;

 j S y , ( ) is the step change indicator for ZPET*;

AD j ; ,S is the step change coefficient for E y 

j S, ( ) for 

the dry period D, principal component j, 

and season S;

AW j ; ,S is the step change coefficient for E y 

j S, ( ) for 

the wet period W, principal component j, 

and season S;

1970 y 2011 I ≤ ≤ is 1 for years 1912–69 and zero otherwise;

1970 y 2011 I ≤ ≤ is 1 for years 1970–2011 and zero otherwise;

P y 

j S, ( ) is the adjusted value of ZP* for the principal 

component j, season S, and year y;

 j S y , ( ) is the step change indicator for ZP*;

BD j ; ,S is the step change coefficient for  j S y , ( ) for 

the dry period D, principal component j, 

and season S; and

BW j ; ,S is the step change coefficient for  j S y , ( ) for 

the wet period W, principal component j, 

and season S.

The values for the step change coefficients (AD;j,s, AW;j,s, BD;j,s, 

BW;j,s) are either zero (if there was no significant change) or 

the mean of the variable for the appropriate time period (if 

the change was significant). Subscript D refers to the “dry” 

period (1912–69) and subscript W refers to the “wet” period 

(1970–2011). A time series model was developed for simulating the adjusted time series. A separate model was used for 

simulating the step change indicators, which were then added 

back into the adjusted time series as described later in this 

section.
Stochastic Climate Model for Simulation of Precipitation, Temperature, and Potential Evapotranspiration 17

For each principal component number (j=1, 2, or 3), the 

six adjusted time series variables (eq. 7) were modeled using a 

multivariate, periodic autoregressive model (McLeod, 1994). 

In this model, the variables were ordered in the following 

sequence:

..., ( ), ( ), ( ), ( ), ( ), , , , , ,

,

    



P y E y P y E y P y

E

j j j j j

j

3 3 1 1 2

2

− − 1 1

( ), ( ), ( ), ( ), ( ),... , , , , y P y E y P y E y j j j j     3 3 1 1 + + 1 1

and each variable was predicted using antecedent values in 

the sequence up to a maximum lag of 1 year (3 seasons); 

for example, P y 

j, ( ) 3 was predicted using potential predictor 

variables P y  E y P y  E y P y  E y j j , , j j , , j j , , ( ), ( ), ( ), ( ), ( ), ( ). 3 3 1 1 2 2 − − 1 1

The coefficients of the prediction equations were assumed to

change from season to season but repeat themselves from year 

to year. The fitted model coefficients are given in table 3. For 

each model equation, a stepwise procedure was used, in which 

the single best antecedent predictor variable was selected. If 

the coefficient for that predictor variable was not significant 

(p≥0.05, where p is the p-value for the coefficient), there were 

no predictor variables included in that model equation, which 

was the case for one equation for PET (PC3, season 1) and 

six equations for precipitation (PC1, season 2; PC2, seasons 1 

and 2; and PC3, all three seasons). The subscripted Z variables 

in the equations are independent, standard normal-random 

variables. If the single best predictor variable was significant 

(p<0.05), the second most significant variable (given the first 

was already included) was selected. If the second variable 

was not significant, the model with the single predictor variable was used. This was the case for five of the nine model 

equations for PET and three of the nine model equations for 

precipitation; for example, for PC1 and season 1, PET was 

negatively related to concurrent precipitation, and precipitation was positively related to precipitation from season 3 of 

the previous year. There were only three equations in which 

two predictor variables were included, all for PET. There were 

no equations with more than two predictor variables. 

Table 3. Fitted multivariate periodic autoregressive models for adjusted principal component time series variables. 

[PCj, jth principal component; E

j S, , adjusted time series variable for potential evapotranspiration for principal component j, season S; Pj S, , adjusted time series 

variable for precipitation for principal component j, season S; y, year; ZE j ; ,S , independent standard normal random variable for potential evapotranspiration principal component j, season S; ZP j ; ,S , independent standard normal random variable for precipitation principal component j, season S]

Principal 

component Season Potential Evapotranspiration Precipitation

PC1 1 E y  

1 1 P y 1 1 ZE 1 1 0 63 1 62 , , ; , ( ) = − . ( ) . + P y  

1 1 P y 1 3 ZP 1 1 0 25 1 1 60 , , ; , ( ) = − . ( ) . +

2 E y  

1 2 P y 1 2 E y 1 1 ZE 1 2 0 30 0 17 1 78 , , , ; , ( ) = − . ( ) . + + ( ) . P y 

1 2 P 1 2 1 71 , ; , ( ) = . Z

3 E y  

1 3 P y 1 3 ZE 1 1 0 41 1 82 , , ; , ( ) = − . ( ) . + P y  P y 1 3 1 2 P 1 3 0 21 1 59 , , ; , ( ) = + . ( ) . Z

PC2 1 E y  

2 1 P y 2 1 ZE 2 1 0 11 0 25 , , ; , ( ) = + . ( ) . P y 

2 1 ZP 2 1 0 75 , ; , ( ) = .

2 E y  

2 2 P y 2 2 E y 2 3 ZE 1 2 0 12 0 19 1 0 28 , , , ; , ( ) = − . ( ) . + − ( ) + . P y 

2 2 ZP 2 2 0 77 , ; , ( ) = .

3 E y  

2 3 E y 2 2 ZE 2 3 0 45 0 33 , , ; , ( ) = + . ( ) . P y  P y 2 3 2 2 P 2 3 0 23 0 79 , , ; , ( ) = + . ( ) . Z

PC3 1 E y 

3 1 ZE 3 1 0 20 , ; , ( ) = . P y 

3 1 ZP 3 1 0 66 , ; , ( ) = .

2 E y  

3 2 P y 3 2 E y 3 3 ZE 3 2 0 13 0 31 1 0 20 , , , ; , ( ) = − . ( ) . + − ( ) + . P y 

3 2 ZP 3 2 0 57 , ; , ( ) = .

3 E y  

3 3 E y 3 2 ZE 3 3 0 67 0 24 , , ; , ( ) = + . ( ) . P y 

3 3 ZP 3 3 0 58 , ; , ( ) = .

PC4 1 E y 

4 1 ZE 4 1 0 09 , ; , ( ) = . P y 4 1 ZP 4 1 0 50 , ; , ( ) = .

2 E y 

4 2 ZE 4 2 0 14 , ; , ( ) = . P y 

4 2 ZP 4 2 0 35 , ; , ( ) = .

3 E y 

4 3 ZE 4 3 0 19 , ; , ( ) = . P y 

4 3 ZP 4 3 0 46 , ; , ( ) = .
18 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

Stochastic Simulation Results

To compare the short-term (annual to decadal) statistical 

properties of the simulated data to the observed data, the time 

series model was used to simulate 100 independent realizations of the time series variables for the dry (1912–69) and wet 

(1970–2011) periods. Each realization was the same length as 

the observed record for the respective dry or wet period. The 

time series model for the principal components (table 3) was 

used to generate the principal component time series. Then, 

the step changes (eq. 7) were added back in for either the dry 

or wet period. The principal components (eq. 3, table 2) were 

reversed (the coefficient matrix for each season was inverted 

and multiplied times the vector of principal components), and 

finally the transformations (eq. 2, table 1) were reversed to 

obtain simulated time series of PET and precipitation for each 

of the station groups and seasons, assuming either dry or wet 

conditions. Statistics computed from the 100 sets of simulated 

data were compared to corresponding statistics of the observed 

record (tables 4 and 5). The statistics included the mean, 10th 

percentile (Q10), and 90th percentile (Q90). The moisture 

deficit (precipitation minus PET) was used in place of PET 

because it is the most relevant variable from a water-balance 

perspective, and it includes the effects of cross-correlation 

between precipitation and PET. The annual total of the three 

seasons also is shown in tables 4 and 5. An 80-percent range 

is shown for the simulated realizations, which is the range 

between the 10th and 90th percentiles out of 100 simulated 

values. 

For most of the statistics, the observed values were well 

within the 80-percent ranges from the simulated realizations, 

and the few that were outside of the ranges were close to the 

upper or lower limits (tables 4 and 5); therefore, the stochastic 

climate model was effective for reproducing the short-term 

variability of the precipitation and PET data for the three seasons and four station groups. Comparison of the dry and wet 

periods provides insight into the potential effects that changes 

in climatic conditions may have on the water-balance model 

simulations described in the “Stochastic Natural Streamflow 

Model” section of this report. For season 1 (November–February), the ranges from the simulations were nearly identical 

between the wet and dry periods. For season 2 (March–June) 

the precipitation statistics tended to be about 10 percent higher 

for all four station groups for the wet period compared to the 

dry period; however, when considering the combined effects 

of precipitation and PET, the deficits were similar for the wet 

and dry periods except for the northeast station group, for 

which deficits tended to be smaller (and surpluses larger) for 

the wet period. For season 3 (July–October), similar to season 

2, the precipitation statistics tended to be about 10 percent 

higher for the wet period; however, unlike season 2, when 

considering the combined effects of precipitation and PET, the 

deficits tended to be lower for the wet period for all four station groups compared to the dry period. These lower summer 

and early fall deficits would be expected to increase runoff 

the following spring as less winter and spring precipitation 

would be required to replace the deficits. Taken in aggregate 

for the annual totals and the deficits, the dry years (indicated 

by Q10) tended to be less severe for the northeast group and 

slightly less severe for the remaining groups for the wet period 

compared to the dry period. The wet years (indicated by Q90) 

tended to be substantially wetter for the southeast, southwest, 

and northeast station groups for the wet period compared to 

the dry period. 

In addition to reproducing the short-term variability 

of the observed climatic data, it is also important that the 

stochastic climate model is able to reproduce the long-term 

(multidecadal to century scale) variability. A 100-year observation period (1912–2011) may be much too short to evaluate 

probabilities of extreme events such as the 2011 spring rainfall 

amounts in the Souris River Basin, especially if climatic conditions are subject to long-range dependence.

Long-term precipitation variability was evaluated using 

12-year moving averages, as was done in the previous section 

of this report “Analysis of Long-Term Climate Variability.” 

An example of simulated 12-year moving average precipitation from the stochastic climate simulation model is shown 

in figure 12 for season 2 (March–June) and station group 1 

(southeast). Although this station group is not the same as 

cluster 1 of the previous analysis (fig. 6), the long-term variability is expected to be similar to cluster 1. Assuming dry 

(1912–69) conditions, for the simulation shown in figure 12A

there were 10 intervals during the 900-year simulation period 

for which the simulated 12-year moving average precipitation 

was less than the recorded low, including three very extreme 

(in both magnitude and duration) droughts comparable to the 

extreme early 1700s drought in figure 6. Conversely, there 

were only two short intervals for which simulated precipitation was higher than the recorded maximum and an 800-

year period during which precipitation remained below the 

recorded maximum. Compared with figure 6, in which there 

were four extreme wet periods (two with long duration) during 

the 300-year extended period, the simulations produce far too 

few extreme wet periods; however, assuming wet (1970–2011) 

conditions (fig. 12B) this simulation produced six wet intervals 

more severe than the observed record and only three brief dry 

intervals comparable to the observed record. 

Comparing figure 6 with figures 12A and 12B, it seems 

that the simulations assuming dry conditions reasonably 

represent the low extremes of the extended record and the 

simulations assuming wet conditions reasonably represent 

the high extremes. Moreover, the extended record was either 

well below or well above the long-term average (and not near 

the long-term average) for most of the 300-year extended 

record; thus, the extended record seems to suggest a mixture 

of alternating dry and wet conditions. As discussed previously 

in this section, strong evidence from the Devils Lake Basin 

in northeastern North Dakota indicates that alternating dry/

wet conditions have been seen in that basin for thousands of 

years (Vecchia, 2008). It is reasonable to postulate that the 

same regional phenomenon is also responsible for alternating 

dry/wet conditions in the Souris River Basin; therefore, the 
Stochastic Climate Model for Simulation of Precipitation, Temperature, and Potential Evapotranspiration 19

Table 4. Statistics of simulated and observed precipitation and moisture deficit (precipitation minus potential evapotranspiration) for 

dry (1912–69) period.

[All units are in millimeters; SE, southeast; SW, southwest; NW, northwest; NE, northeast; Annual, total of all three seasons; top number is for the observed 

record; numbers below represent the 10th and 90th percentiles from 100 simulated datasets]

Season Group Precipitation Precipitation minus potential 

evapotranspiration

Mean 10th percentile 90th percentile Mean 10th percentile 90th percentile

1

SE 53 32 86 18 -8 56

46, 52 24, 31 66, 80 10, 18 -17, -6 33, 50

SW 50 31 69 13 -9 40

48, 53 27, 34 66, 80 9, 17 -17, -8 32, 46

NW 98 65 131 69 34 104

89, 97 58, 68 116, 132 60, 68 24, 35 88, 108

NE 84 58 119 53 18 91

78, 87 50, 58 106, 123 46, 56 13, 24 77, 94

2

SE 180 107 255 -31 -112 52

171, 188 104, 127 229, 263 -42, -22 -116, -89 23, 61

SW 167 98 230 -50 -132 29

157, 175 90, 116 215, 252 -62, -41 -136, -107 2, 42

NW 183 115 264 -9 -80 69

170, 194 101, 128 238, 277 -23, 2 -102, -70 50, 93

NE 195 122 268 -14 -91 47

187, 206 117, 143 247, 283 -24, -1 -100, -71 40, 82

3

SE 180 125 234 -122 -190 -51

173, 192 112, 132 228, 257 -132, -110 -205, -177 -68, -36

SW 148 101 186 -159 -223 -98

141, 154 90, 107 186, 212 -167, -153 -230, -211 -112, -86

NW 183 113 262 -89 -167 13

173, 187 110, 132 227, 258 -102, -85 -172, -147 -37, -5

NE 206 143 254 -93 -157 -33

195, 215 132, 156 255, 289 -107, -85 -182, -153 -35, -1

Annual

SE 413 315 495 -135 -227 -32

392, 424 287, 332 484, 534 -160, -123 -274, -231 -53, 2

SW 365 273 451 -196 -296 -91

351, 377 258, 294 434, 478 -216, -183 -320, -281 -115, -67

NW 465 348 597 -28 -150 111

441, 473 331, 374 543, 596 -53, -19 -178, -130 59, 123

NE 485 364 599 -54 -178 84

468, 500 356, 405 569, 617 -75, -40 -197, -149 35, 91
20 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

Table 5. Statistics of simulated and observed precipitation and moisture deficit (precipitation minus potential evapotranspiration) for 

wet (1970−2011) period.

[All units are in millimeters; SE, southeast; SW, southwest; NW, northwest; NE, northeast; Annual, total of all three seasons; top number is for the observed 

record; numbers below represent the 10th and 90th percentiles from 100 simulated datasets]

Season Group Precipitation Precipitation minus potential 

evapotranspiration

Mean 10th percentile 90th percentile Mean 10th percentile 90th percentile

1

SE 49 22 75 12 -16 34

46, 53 24, 32 67, 83 8, 16 -19, -8 33, 49

SW 56 33 87 16 -14 53

47, 54 27, 33 66, 79 6, 14 -21, -9 28, 44

NW 97 65 133 65 26 103

89, 98 58, 70 116, 134 57, 68 20, 36 88, 106

NE 91 59 134 57 20 105

79, 88 49, 60 104, 127 44, 54 9, 22 72, 96

2

SE 193 127 271 -25 -111 54

184, 206 114, 143 242, 288 -35, -11 -113, -83 30, 82

SW 178 115 251 -49 -118 32

170, 193 102, 130 228, 272 -58, -33 -132, -103 3, 56

NW 200 138 274 -6 -74 71

189, 213 117, 144 251, 295 -19, 8 -99, -67 54, 102

NE 216 164 301 -3 -62 86

202, 225 129, 158 265, 308 -16, 9 -96, -66 52, 101

3

SE 204 135 268 -96 -172 -28

190, 212 125, 154 250, 286 -110, -86 -185, -153 -42, -1

SW 164 108 226 -148 -228 -79

154, 172 102, 126 201, 230 -158, -137 -224, -197 -103, -68

NW 196 143 254 -82 -141 -24

186, 208 124, 153 243, 282 -91, -67 -164, -128 -32, 13

NE 227 139 299 -71 -176 10

216, 240 150, 178 277, 316 -81, -56 -162, -125 -12, 35

Annual

SE 447 348 549 -110 -226 21

430, 464 316, 365 526, 581 -129, -90 -249, -199 -23, 39

SW 398 302 508 -181 -279 -63

377, 414 282, 325 460, 518 -203, -164 -316, -264 -100, -36

NW 493 397 615 -23 -133 85

473, 510 361, 411 564, 643 -43, -2 -169, -107 60, 142

NE 533 436 647 -17 -127 73

508, 542 395, 441 605, 674 -43, -4 -174, -121 69, 131
Stochastic Climate Model for Simulation of Precipitation, Temperature, and Potential Evapotranspiration 21

A. Dry (1912–69) climatic conditions

120

140

160

180

200

220

240

B. Wet (1970–2011) climatic conditions

120

140

160

180

200

220

240

Precipitation, in millimeters

C. Alternating dry (unshaded) and wet (shaded) climatic conditions

120

140

160

180

200

220

240

0 100 200 300 400 500 600 700 800 900

Simulation year

Maximum observed value, 1912–2011

Minimum observed value, 1912–2011

Maximum observed value, 1912–2011

Minimum observed value, 1912–2011

Maximum observed value, 

1912–2011

Minimum observed 

value, 1912–2011

Figure 12. Long-term (900-year) simulated 12-year moving average precipitation for season 2 (March–June) and station 

group 1 (southeast). A, dry conditions. B, wet conditions. C, alternating dry/wet conditions.
22 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

same alternating long-term simulated climate states that were 

used for the Devils Lake simulation model (Vecchia, 2008), 

in which dry periods averaged 120 years in duration and wet 

periods 30 years in duration, were used to model the shifting 

climate states for the Souris River Basin. An example of simulated 12-year moving average precipitation with alternating 

dry/wet climate states is shown in figure 12C. In this example, 

there were six wet states with durations ranging from about 

10 to 60 years interspersed with dry states ranging from about 

50 to more than 200 years in duration. There were 9 severe 

droughts, 8 of which took place during dry states, and 6 major 

wet intervals (above 220 millimeters [mm]), 5 of which took 

place during wet states.

An example of simulated 12-year moving average 

precipitation from the stochastic climate simulation model 

is shown in figure 13 for season 3 (July–October) and station group 4 (northeast). Although this station group is not 

the same as cluster 5 of the previous analysis (fig. 7), the 

long-term variability is expected to be similar to cluster 5. 

Similar to the previous example with station group 1 (southeast) and season 2 (March – June), the simulation assuming 

dry conditions (fig. 13A) reproduces the extreme droughts of 

the extended record (fig. 7), and the simulation assuming wet 

conditions (fig. 13B) reproduces the extreme wet periods. The 

simulation with alternating dry/wet conditions (fig. 13C) has 9 

severe droughts, 8 of which took place during dry states, and 4 

prolonged wet periods (above 240 mm), 3 of which took place 

during wet states.

Water-Balance Model for Estimating 

Natural Streamflow

To convert stochastic estimates of precipitation, temperature, and PET into streamflow for flood risk analysis, 

evaluation of reservoir operation, and general simulation 

of future streamflow, a water-balance model (WBM) was 

developed such that three variables (precipitation, temperature, and PET), in conjunction with static basin properties (soil 

water storage and permeability), could be used to estimate 

monthly streamflow. The model was calibrated with historic 

(1946–2011) average monthly natural (unregulated) streamflow data provided by the U.S. Army Corps of Engineers 

(USACE) and with observed (gaged) streamflow for locations 

where observed flow was the same as natural flow. The natural 

streamflow data provided by the USACE was unregulated, 

which means the effects of Alameda, Rafferty, Boundary, and 

Lake Darling reservoirs/dams were removed; and streamflow 

was allowed to flow through the basin as it would have before 

the four reservoirs were installed. Once the WBM had been 

calibrated with the monthly natural streamflow, flows were 

then disaggregated into an approximately 10-day time step 

(three values per month). 

Input Data for Water-Balance Model

The WBM requires three dynamic inputs (monthly 

precipitation, average temperature, and PET) and two static 

inputs (soil water storage and permeability) to simulate natural 

streamflow at selected streamflow-gaging stations (hereafter 

referred to as “gaging stations”). In order to determine the 

contributing area to selected gaging stations, basins were 

delineated using pre-defined basin delineations, online delineation tools, and (or) digital elevation models (DEMs). Calibration of the model relied heavily on streamflow estimates 

at various gaging stations throughout the Souris River Basin; 

therefore, selection of gaging stations with long and continuous streamflow records was critical.

Streamflow

Gaging stations were selected based on location and 

length of historical streamflow record. Streamflow data were 

retrieved from U.S. Geological Survey (2014b) and Water 

Survey of Canada (Government of Canada, 2014) online databases. A list of all gaging stations used to develop the WBM is 

provided in table 6 and figure 1.

In addition to observed streamflow, natural streamflow 

estimates from the USACE (U.S. Army Corps of Engineers, 

2013) were used to calibrate the WBM. Natural streamflow 

is that streamflow that would have naturally occurred without the installation of the Boundary, Rafferty, Alameda, and 

Lake Darling Reservoirs (fig. 1). Natural streamflow does 

take into account the effects of the Des Lacs and J. Clark 

Salyer National Wildlife Refuges (NWRs; fig. 1), along with 

the effects of other minor hydraulic structures throughout the 

Souris River Basin, because these structures allow water to 

pass in a way that mimics natural streamflow patterns and 

are not operated for flood mitigation (U.S. Army Corps of 

Engineers, 2013). Natural flows were estimated through the 

use of five different flow routing methods: Null routing (direct 

routing), Muskingum routing, Muskingum Cunge routing, 

Lag routing, and Modified Puls routing (U.S. Army Corps of 

Engineers, 2013). 

Watershed Delineation

Selected gaging stations were designated as outlet points, 

and contributing drainage areas were delineated resulting in 

the formation of 19 subbasins (both intervening and headwater basins) within the Souris River Basin (fig. 14). The four 

Manitoba gaging stations (05NF001, 05NG021, 05NG007, 

and 05NG001; fig. 14) downstream from the United States are 

not included in the following model development, but basin 

delineations for these gaging stations were developed for 

potential future extensions of the model from the international 

border to the basin outlet.
Water-Balance Model for Estimating Natural Streamflow 23

A. Dry (1912–69) climatic conditions

260

280

160

180

200

220

240

260

280

160

180

200

220

240

260

280

160

180

200

220

240

B. Wet (1970–2011) climatic conditions

Precipitation, in millimeters

C. Alternating dry (unshaded) and wet (shaded) climatic conditions

0 100 200 300 400 500 600 700 800 900

Simulation year

Maximum observed value, 1912–2011

Minimum observed 

value, 1912–2011

Maximum observed 

value, 1912–2011

Maximum observed 

value, 1912–2011

Minimum observed value, 1912–2011

Minimum observed 

value, 1912–2011

Figure 13. Long-term (900-year) simulated 12-year moving average precipitation for season 3 (July–October) and station 

group 4 (northeast). A, dry conditions. B, wet conditions. C, alternating dry/wet conditions.
24 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

Table 6. Selected United States and Canadian streamflow-gaging stations used for calibration of the water-balance model.

[Sask., Saskatchewan; N. Dak., North Dakota; Man., Manitoba]

Station

identifier Station name State/province

Period of 

record Comment

05ND004 Moose Mountain Creek near Oxbow Sask. 1913−present Seasonal with some continuous 

record. No record 1918−32.

05NA003 Long Creek at western crossing of international 

boundary

Sask. 1959−present Continuous.

05113600 Long Creek near Noonan N. Dak. 1959−present Continuous.

05NB001 Long Creek near Estevan Sask. 1911−present Seasonal and continuous. No record 

1924−32 and 1958.

05NB036 Souris River below Rafferty Reservoir Sask. 1992−present Continuous.

05114000 Souris River near Sherwood N. Dak. 1930−present Continuous.

05116000 Souris River near Foxholm N. Dak. 1936−present Continuous.

05116500 Des Lacs River at Foxholm N. Dak. 1904−present No record 1907−45.

05117500 Souris River above Minot N. Dak. 1903−present Continuous.

05120000 Souris River near Verendrye N. Dak. 1937−present Continuous.

05120500 Wintering River near Karlsruhe N. Dak. 1937−present Continuous.

05122000 Souris River near Bantry N. Dak. 1937−present Continuous.

05123400 Willow Creek near Willow City N. Dak. 1956−present Continuous.

05123510 Deep River near Upham N. Dak. 1957−present Seasonal.

05124000 Souris River near Westhope N. Dak. 1929−present Continuous.

05NF001 Souris River at Melita Man. 1912−present Seasonal with some continuous 

record. No record 1923−34, 

1937−74, 1977−2001.

05NG007 Plum Creek near Souris Man. 1956−present Seasonal with some continuous record. No record 1995−2001.

05NG021 Souris River at Souris Man. 1946−present Seasonal with some continuous 

record. No record 1948−66.

05NG001 Souris River at Wawanesa Man. 1912−present Seasonal with some continuous 

record.

An order of operations was followed when determining 

the method used for subbasin delineation. First, the Souris 

River Basin delineation, provided through the Agriculture 

and Agri-Food Canada’s (AAFC) Watersheds Project 2013 

(Government of Canada, 2015b), defined the border/edge of 

the basin for the Souris River. Second, USGS StreamStats, 

an online watershed delineation tool (U.S. Geological Survey, 2014c), was used for delineation of subbasins with an 

outlet located in the United States. To check the accuracy of 

the StreamStats delineation, and assist with determining the 

contributing area, the USGS Watershed Boundary Dataset 

(U.S. Geological Survey, 2015) was used. For instances when 

gaps between StreamStats and the USGS Watershed Boundary 

Dataset existed, delineation with ArcMap 10.2.2 Hydro Tools 

(Esri, 2014) and USGS 10-meter DEM (U.S. Geological Survey, 2014a) data was used. For outlets in Canada, the AAFC 

Total Gross Drain dataset (Government of Canada, 2015a) was 

used to define catchment areas given selected gaging stations. 

Soil Characteristics

Soil permeability and available water storage data were 

obtained from the Gridded Soil Survey Geographic Database 

(gSSURGO) version 1.1 (U.S. Department of Agriculture, 

2014), Detailed Soil Survey (DSS) of Canada (Agriculture 

and Agri-Food Canada, 2014), and Soil Landscapes Canada 

(SLC) Version 3.2 (Agriculture and Agri-Food Canada, 2013) 

at scales of 1:24,000; 1:100,000; and 1:1 million, respectively. 

Weighted averages were used to estimate soil permeability and 

available water storage capacity over 90- and 100-centimeter 

depths, respectively. 

Weather Data

Monthly precipitation and temperature data were 

retrieved from the USHCN and AHCCD online databases as 

described in the previous section “Analysis of Long-Term 

Climate Variability.” A total of 27 sites were used to estimate 
Water-Balance Model for Estimating Natural Streamflow 25

! !

!! !

!

!!

!

! !

!

!

!

!

!

!! !

!

Ç

Ç

Ç

Ç

Ç

Ç

Ç Ç

Ç

Ç

Ç

Ç

Ç

Ç

Ç

Oxbow

Minot

Upham

Souris

Me lita

Noona n

Bantry

Estevan

Foxholm

Wawa ne sa

She rwood

We sthope

Verendrye

Karlsruhe

Willow City

Ç

Ç

Ç

Minot

Verendrye

Karlsruhe

( (

(( (

(

((

(

( (

(

(

(

(

(

(( (

99° 100° 101° 102° 103° 104° 105°

50°49°

48°

(

MANITOBA

SASKATCHEWAN

MONTANA

NORTH DAKOTA

CANADA

UNITED STATES

Intervening basin (flow routing required)Headwater basin (no flow routing required)

Aquifer (index map)

EXPLANATION

KarlsruheNew RockfordMinot

Streamflow-gaging station

SundreVoltaire

0 60 80 120 KILOMETERS 300 40 60 MILES 20

 Index map

Figure 14.balance model.

Streamflow-gaging stations and corresponding subbasins, designated by station identifier (table 1), used for developing the water-
26 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

monthly precipitation, temperature, and PET throughout the 

basin (table 7), with 26 of those sites located within a 1-degree 

buffer of the Souris River basin (fig. 8). Monthly PET values 

were estimated through use of the Hamon method, as previously described in the section “Stochastic Climate Model 

for Simulation of Precipitation, Temperature, and Potential 

Evapotranspiration.”

In comparison with other PET estimation methods, the 

Hamon method has been shown to underestimate monthly 

PET, especially during spring and early summer (Rao and others, 2011). Rao and others (2011) suggested underprediction of 

PET by the Hamon method may be attributed to the method’s 

lack of consideration of solar radiation. To compensate for 

low PET values in the WBM, a correction factor was applied 

such that PET values for every month were increased by the 

same percent. Through trial and error, a 10-percent increase 

in monthly PET was determined to be effective for the WBM. 

This correction factor seemed reasonable given other applications, such as Lu and others (2005), which used a correction 

factor to increase monthly PET by 20 percent. Similarly, a 

study by Rosenberry and others (2004) determined slight 

increases in PET, estimated using the Hamon method, resulted 

in better comparisons with the Bowen-ratio energy-budget 

method they used as standard in the study. In addition to the 

general underprediction of PET using the Hamon method, Rao 

and others (2011) noted the Hamon method had a tendency to 

stray from seasonal evapotranspiration (ET) patterns compared 

to other methods (such as Priestley-Taylor and Food and Agriculture Organization [FAO] full Penman-Monteith equation 

for a grass surface), with PET estimates that were significantly 

lower during spring and early summer. This low PET in the 

spring and early summer led to a steeper ascent to peak PET 

in July. To compensate for this underprediction and steep 

ascent, PET values from the Hamon method were increased 

in May and June, such that each subbasin adjustment could 

be different. By increasing the PET values, the steep ascent to 

peak PET in July was reduced along with reducing the general 

underprediction of PET. A list of watershed-specific PET 

adjustments is provided in table 8.

 Water-Balance Model Description and 

Calibration 

The WBM, adapted from Gray and McCabe (2010), is a 

mass-based deterministic model. Climatic inputs to the model 

consist of monthly precipitation, average temperature, and 

PET, along with static inputs such as soil permeability and 

available water storage capacity measured in depth equivalents. Climatic inputs are routed through the model using a 

series of equations intended to estimate the amount of water 

added, retained, and (or) lost at each monthly time step (see 

appendix). Application of the WBM at the scale of the subbasins (fig. 14) would have limited the ability of the model 

to take into account changes in soils and climatic variability 

within the subbasins. To better account for variations in soils 

and climatic conditions throughout the subbasins, and at the 

same time maintain reasonable computation times, a grid 

of 8x8-km cells was used to divide the Souris River Basin 

into 1,103 cells. The WBM was then applied to individual 

grid cells, allowing for site-specific estimation of runoff in 

response to spatial variations in precipitation and (or) soil 

characteristics across the basin. Soil permeability and available water storage capacity were averaged to each grid cell 

using weighted averages (fig. 15). Precipitation, temperature, 

and PET were interpolated from the weather station locations 

(fig. 8) to each grid cell using locally weighted polynomial 

regression (loess; Cleveland and others, 1992) in the statistical 

Table 7. Selected United States and Canadian weather stations 

used for calibration of the water-balance model.

[Sask., Saskatchewan; Mont., Montana; N. Dak., North Dakota; Man., Manitoba]

Station

identifier Station name State/province Latitude Longitude

4011441 Ceylon Sask. 49.4 -104.7

4016560 Regina Sask. 50.4 -104.7

245572 Medicine Lake 3 

southeast

Mont. 48.5 -104.5

4019040 Yellow Grass Sask. 49.8 -104.2

4013660 Kelliher Sask. 51.3 -103.8

4013480 Indian Head Sask. 50.5 -103.7

321871 Crosby N. Dak. 48.9 -103.3

4012400 Estevan Sask. 49.2 -103.0

4010400 Bangor Sask. 50.9 -102.3

4019082 Tonkin Sask. 51.2 -102.2

4014913 Manor Sask. 49.6 -102.1

4015360 Moosomin Sask. 50.1 -101.7

5012080 Pierson Man. 49.2 -101.3

5010240 Birtle Man. 50.4 -101.1

325479 Mandan Experiment 

Station

N. Dak. 46.8 -100.9

320941 Bottineau N. Dak. 48.8 -100.4

328792 Towner 2 northeast N. Dak. 48.4 -100.4

329445 Willow City N. Dak. 48.6 -100.3

5010480 Brandon Man. 49.9 -99.9

5022000 Ninette Man. 49.4 -99.7

5042004 Neepawa Murray Man. 50.2 -99.6

5010640 Cypress River Man. 49.5 -99.1

5041535 Langruth West Man. 50.4 -98.8

324958 Langdon Experimental FarmN. Dak. 48.8 -98.3

5012321 Portage Prairie Man. 50.0 -98.3

5021849 Morden Man. 49.2 -98.1

248569 Vida Mont. 47.9 -105.4
Water-Balance Model for Estimating Natural Streamflow 27

software R (R Core Team, 2015). The loess function makes 

use of locally weighted regression to assign values to each 

grid cell based on latitude and longitude. For this application, 

similar to Ryberg and others (2014), the loess function with 

polynomial degree of 1 (linear in latitude and longitude) and a 

span of 0.7 was used. 

A simplified schematic of the WBM, adapted from 

Gray and McCabe (2010), for a generic grid cell is shown in 

figure 16. A brief description of the model is given here and 

detailed equations are given in the appendix. Precipitation for 

a grid cell is partitioned into rain or snow depending on the 

average monthly temperature. Snow is stored in the snowpack 

until early spring when temperatures rise and snowmelt is 

produced. A small portion of the snowmelt directly enters the 

stream as runoff. The remainder of the snowmelt is used to 

help reduce (or eliminate) the soil moisture deficit remaining 

from the previous time step and, depending on time of year, 

contributes to subsurface flow. The rate of subsurface flow (the 

percentage of the soil moisture that is assumed to enter the 

stream) increases with increasing permeability and decreases 

with increasing soil moisture deficit. The rate of subsurface 

flow is also dependent on temperature such that temperatures below -10 degrees Celsius (o

C) restrict all subsurface 

flow because soils are frozen. Temperatures above 2 oC leave 

subsurface flow unrestricted because soils are thawing or have 

thawed. After losses to subsurface flow have been considered 

and additions either through snowmelt and (or) liquid precipitation have been accumulated, a fraction of the water in excess 

of the soil moisture capacity (SMC) is considered to contribute 

to direct runoff (if any) and is therefore unavailable for PET. 

The remaining water in the system, after subtractions from 

direct runoff, is referred to as the total available water (TAW) 

(or the water available for PET or additional overland flow). 

If TAW exceeds SMC plus moisture demand (PET) for the 

month, ET is assumed to equal PET, and the excess (TAW–

PET–SMC) is assumed to enter the stream as overland flow. 

The proportion of overland flow entering the stream in the current month, compared to the following month, was determined 

Table 8. Parameters from calibration of the monthly water-balance model using reconstructed natural (unregulated) monthly 

streamflow from the U.S. Army Corps of Engineers.

[CAWS, available water storage capacity parameter; CDRO, direct runoff parameter; CSM, snowmelt runoff parameter; PET, potential evapotranspiration; m3/s, cubic 

meter per second; Sask., Saskatchewan; N. Dak., North Dakota; NA, correlation coefficient not available due to lack of reconstructed streamflow data; --, no 

runoff cap applied; calibration coefficients for each subbasin are described in the appendix]

Station/

subbasin 

identifier

Station name State/

province CAWS CDRO CSM

Adjustment 

factors Correlation 

coefficient1

Runoff 

cap

(m3

/s)

Percentage of 

excess overland 

flow entering 

the stream in the 

current/following 

month

PET 

(May)

PET 

(June)

05ND004 Moose Mountain Creek 

near Oxbow

Sask. 1.2 0.15 0.01 1.35 1.25 0.85 56 50/50

05NA003 Long Creek at western 

crossing of international 

boundary

Sask. 1.0 0.6 0.05 1.2 1.1 0.77 -- 50/50

05113600 Long Creek near Noonan N. Dak. 1.0 0.2 0.03 1.2 1.1 0.78 -- 50/50

05NB001 Long Creek near Estevan Sask. 1.0 0.2 0.03 1.2 1.1 NA -- 50/50

05NB036 Souris River below Rafferty 

Reservoir

Sask. 1.1 0.35 0.01 1.2 1.1 0.78 -- 50/50

05114000 Souris River near Sherwood N. Dak. 1.2 0.1 0.01 1.3 1.2 0.85 -- 50/50

05116000 Souris River near Foxholm N. Dak. 1.0 0.3 0.03 1.3 1.2 0.85 -- 50/50

05116500 Des Lacs River at Foxholm N. Dak. 1.0 0.3 0.04 1.0 1.0 0.75 -- 50/50

05117500 Souris River above Minot N. Dak. 1.0 0.3 0.04 1.0 1.0 0.84 -- 50/50

05120000 Souris River near Verendrye N. Dak. 0.9 0.5 0.05 1.0 1.0 0.84 -- 50/50

05120500 Wintering River near 

Karlsruhe

N. Dak. 1.0 0.5 0.04 1.1 1.0 0.75 21 50/50

05122000 Souris River near Bantry N. Dak. 1.0 0.3 0.03 1.1 1.0 0.84 -- 50/50

05123400 Willow Creek near Willow 

City

N. Dak. 1.0 0.25 0.01 1.2 1.1 0.81 -- 50/50

05123510 Deep River near Upham N. Dak. 1.1 0.4 0.01 1.1 1.0 0.83 -- 80/20

05124000 Souris River near Westhope N. Dak. 1.1 0.1 0.01 1.1 1.0 0.78 -- 20/80

1

 Correlation coefficient is of log transformed data. Before data was logged each streamflow value was increased by 1 m3/s to eliminate the possibility of taking the log of a zero flow.
28 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

±

0 –70 71–90 91–110 111–130 131–150 151–170 171–190 191– 210 211– 230 231– 250

0– 2 3 – 4 5– 6 7– 8 9 –10 11–12 13 –14 15–16 17 –18 19 –20 21+

Available water storage, in millimeters

Soil permeability, in centimeters per hour

0 60 80 120 KILOMETERS 300 40 60 MILES 20

CANADA

UNITED STATES

CANADA

UNITED STATES

CANADA

UNITED STATES

CANADA

UNITED STATES

A

C

BD

Figure 15. Soil permeability and available water storage (AWS) capacity. A, soil permeability before data are averaged to the 8-kilometer x 8-kilometer grid. 

B, soil 

permeability after data are averaged to the 8-kilometer x 8-kilometer grid. 

C, AWS before data are averaged to the 8-kilometer x 8-kilometer grid. 

D, AWS after data are 

averaged to the 8-kilometer x 8-kilometer grid.
Water-Balance Model for Estimating Natural Streamflow 29

Precipitation Runoff

Temperature

Permeability

Primary ET

min(SMS, PET)

ET PET

Secondary ET

~ SM, SMC

Snow Snowpack

Rain

Snowmelt

Soil moisture

surplus (SMS)

Soil moisture (SM)

Soil moisture capacity (SMC)

Early

snowmelt

runoff

Overland

flow

Excess

snowmelt

runoff

Subsurface

flow

Evapotranspiration

EXPLANATION

Potential evapotranspiration

min() Minimum

Figure 16. Schematic of water-balance model for a generic grid 

cell.

through trial and error (table 8) with the intent to closely 

match the monthly means and standard deviations of the 

natural streamflow reported by the USACE (U.S. Army Corps 

of Engineers, 2013). Basins in other settings would most likely 

have different proportions of flow entering the stream in the 

current as compared to the following month, because the proportion is highly dependent on size, topography, and climate 

of the basin. In addition, for 2 (Moose Mountain Creek near 

Oxbow, Sask., 05ND004 and Wintering River near Karlsruhe, 

N. Dak., 05120500) of the 15 subbasins, the amount of overland flow was capped such that the overland flow for a given 

month could not exceed the maximum streamflow recorded 

(1946 – 2011) at the outlet of that basin (table 8). Runoff caps 

were applied to those subbasins whose monthly mean streamflow values, after calibration, still tended to overpredict the 

monthly mean streamflow values provided by the USACE 

(U.S. Army Corps of Engineers, 2013).

If TAW exceeds soil moisture capacity but is not sufficient to satisfy PET, the surplus (TAW–SMC) is lost as 

primary ET, and there is no overland flow. A part of the unsatisfied moisture demand (PET–primary ET) is subtracted from 

the soil moisture as secondary ET, where the amount of secondary ET decreases in proportion to the ratio of soil-moisture 

storage to the soil-moisture storage capacity. If TAW is less 

than SMC, there is no primary ET, and thus total ET consists 

of secondary ET removed from the soil moisture. 

The model has five adjustable parameters. These five 

parameters were estimated to maximize the correlation 

between estimated and natural monthly streamflows while 

matching, as closely as possible, monthly means and standard 

deviations of the natural monthly streamflows (table 8). Two 

of the five parameters (FGW and FST in the appendix) were 

dependent on spatially variable soil properties and temperature 

and were not calibrated to individual subbasins. The three 

basin-specific parameters consist of two parameters to control 

the amount of early snowmelt runoff and the rate of surplus 

runoff and one parameter to adjust the available water storage capacity. The snowmelt runoff parameter (CSM) ranged 

between 0.01 and 0.05, the surplus runoff parameter (CDRO) 

between 0.1 and 0.6, and the available water storage capacity 

parameter (CAWS) between 0.9 and 1.2 (table 8). The parameters act as coefficients within the model, either increasing or 

decreasing model values by a certain percent. Neighboring or 

upstream subbasins tended to have similar parameters as seen 

with the Souris River near Foxholm, N. Dak. (05116000); Des 

Lacs River at Foxholm, N. Dak. (05116500); and Souris River 

above Minot, N. Dak (05117500) subbasins (fig. 14), where 

each basin had a CAWS andCDRO value of 1.0 and 0.3, respectively, and varied slightly with CSM values ranging between 

0.03 and 0.04 (table 8). Subbasins that tended to have higher 

CAWS values also tended to have lower CDRO and CSM values. 

Similarly, subbasins with lower CAWS tended to have higher 

CDRO and CSM values. This can be seen with the Moose Mountain Creek near Oxbow, Sask. (05ND004) subbasin, which had 

CAWS, CDRO, and CSM values of 1.2, 0.15, and 0.01, respectively; 

and the Souris River near Verendrye, N. Dak. (05120000) 

subbasin, which had CAWS, CDRO, and CSM values of 0.9, 0.5, and 

0.05, respectively (table 8). 

The WBM estimates monthly runoff for an individual 

grid cell; taking the average of all grid cells overlying a 

specific subbasin provides an average runoff estimate for that 

specific subbasin. Runoff was estimated for each of the 15 

subbasins; subbasins that did not require streamflow routing (8 

headwater basins) were calibrated first through trial and error 

(fig. 14). Two objectives were used in the calibration: maximizing the correlation between estimated and natural monthly 

flows and reproducing, as closely as possible, the monthly 

means and standard deviations of the natural flows. The 

WBM is not intended as a forecast model to predict natural 

flows as closely as possible. A forecast model would require 

much more detailed data inputs with respect to the spatial and 

temporal distribution of precipitation in the basin and was 

considered beyond the scope of this report. Rather, the WBM 

is intended for simulating flows that reproduce the statistical 

distributions of natural monthly flows given monthly climatic 

inputs. In most cases, a higher correlation between estimated 

and natural flows could have been obtained at the expense of 

deviating from the monthly means and standard deviations; 

furthermore, there may be many different model configurations and (or) parameter estimates that result in similar correlation coefficients. The monthly means and standard deviations were more sensitive to the model equations (appendix) 
30 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

and parameter values (table 8) than were the correlation 

coefficients between estimated and natural flows.

An example of the calibrated WBM output for one of the 

headwater basins, Long Creek at western crossing of the international boundary (05NA003) (fig. 14), is shown in figures 

17 and 18. The WBM was able to closely replicate historical 

streamflow events where simulated results were able to catch 

most of the high flow and almost all of the low flow events 

with a correlation coefficient of 0.8 (fig. 17). In addition, the 

monthly means showed little to no bias between natural and 

simulated values (fig. 18A). Similarly, monthly standard deviations for simulated and natural flows followed the same trend 

with monthly variation peaking in April and June along with a 

slight dip in May (fig. 18B). 

Natural (Unregulated) Streamflow Routing

After runoff values were estimated for each of the 

headwater basins, flows were routed through downstream 

basins and intervening flows calibrated with USACE natural flow data (U.S. Army Corps of Engineers, 2013). During 

the calibration process, intervening basin parameters were 

checked for similarity to neighboring basins to ensure calibration was reasonable (table 8). All but one subbasin, Long 

Creek near Estevan, Sask. (05NB001), was calibrated with 

natural flow data (fig. 14). For the Long Creek near Estevan, 

Sask. (05NB001) subbasin, there was no natural flow data 

available for calibration, and as a result basin parameters were 

assumed to be the same as the upstream subbasin, Long Creek 

near Noonan, N. Dak. (05113600). Given the relatively small 

size of the Long Creek near Estevan, Sask. (05NB001) subbasin (295 km2

) and similar topography and vegetation to the 

Long Creek near Noonan, N. Dak. (05113600) subbasin, this 

seemed to be a reasonable assumption.

To route monthly flows through the Souris River Basin, 

direct routing and lagging were used. Direct routing involved 

the simple addition of streamflow from separate subbasins 

(upstream and downstream) such that flows in the same time 

step were added. Lagging involved delaying the amount of 

flow for the current time step by partitioning the flow such 

that a fraction of the flow comes from the current time step 

and the remaining fraction comes from the previous time step; 

for example, taking 90 percent of the flow from the current 

time step and 10 percent from the previous time step creates a lagged/delayed flow of about 3 days when considering 

a monthly time step of 30 days. The combination of direct 

and lagged routing allows for a simplified routing model that 

takes into consideration time of travel for routed flows. When 

the reach between a basin’s outlet and the downstream outlet 

was short (less than about 45 km), no lag was applied and 

flows were directly routed. Subbasins that did not have a lag 

for downstream flow routing included Long Creek at western crossing of the international boundary (05NA003); Long 

Creek near Noonan, N. Dak. (05113600); Moose Mountain 

Creek near Oxbow, Sask. (05ND004); Souris River near 

Foxholm, N. Dak. (05116000); Des Lacs River at Foxholm. 

N. Dak. (05116500); and Souris River near Westhope, N. Dak. 

(05124000) (table 9). In cases where the river reach extended 

beyond 45 km and (or) there were hydraulic structures limiting flow, lags were applied. The extent of the lag applied 

was determined through trial and error such that the monthly 

means and standard deviations of simulated routed monthly 

flows closely matched those provided by the USACE (U.S. 

Army Corps of Engineers, 2013). A list of the subbasins and 

their respective lags are given in table 9. Lags were applied 

to their respective subbasin when flows were being routed 

from that subbasin to the next (table 9). For example, for 

estimating streamflow at the Souris River near Foxholm, N. 

Dak. (05116000) gaging station, flows (with their respective 

lags) would be routed up to Foxholm, but the lag provided 

for Foxholm would not be applied; however, for estimating 

streamflow at the next downstream gaging station, Souris 

River above Minot, N. Dak. (05117500), the lag would be 

applied to Foxholm as accumulated flows are routed from 

Foxholm to Minot. Simulated monthly means for the Souris 

River near Sherwood, N. Dak. (05114000) and Souris River 

near Westhope, N. Dak. (05124000) gaging stations compared 

well with mean natural monthly streamflow from the USACE 

(fig. 19), showing very little bias between simulated and 

natural monthly streamflow. Similarly, the standard deviations 

for simulated and natural monthly streamflow compared well 

(fig. 20) with simulated monthly standard deviations following 

the same trends as those produced by the USACE data.

Disaggregation of monthly streamflow to a three-permonth or approximately 10-day time step allows for a more 

precise routing of flows and provides simulated natural flow 

data that can assist with reservoir operation. Reservoir operation is highly dependent on 30-day forecasted streamflow volumes, with forecasts being updated daily during flood years. 

Flood years are those in which the 30-day streamflow volume 

at the Souris River near Sherwood, N. Dak. (05114000) gaging station (fig. 1) is greater than 216,105,696 cubic meters 

(m3

) (175,200 acre-foot [acre-ft]) or local flows to this gage 

exceed a 30-day volume of 37,004,400 m3

 (30,000 acre-ft) 

(U.S. Army Corps of Engineers, 2013). Disaggregation of 

monthly streamflow to a three-per-month time step allows 

reservoir operators to not only look at forecasted streamflow volumes for each month but also allows them to look at 

forecasted streamflow volumes for the 1st (days 1–10), 2nd

(days 11–20), and 3rd (days 21 through the end of the month) 

periods of each month, which enables them to make more time 

sensitive operation decisions.

Disaggregation of monthly flows was done by taking 

historic daily flow values at each of the gaging points and 

calculating the fraction of monthly flow occurring in the 1st, 

2nd, and 3rd periods of the month. The disaggregation resulted 

in there being three periods for each month, each of which 

had a ratio defining the fraction of monthly flow occurring 

during that period, such that when the ratios were added for all 

three periods the sum was equal to one. The ratios were then 

used to subdivide simulated monthly streamflow to threeper-month values during the same historic record. In other 
Water-Balance Model for Estimating Natural Streamflow 31

Year

Streamflow, in cubic meters per second

1955 1960 1965 1970 1975 1980 1985 1990 1995 2000 2005 2010 2015

0

10

20

30

40

50

60

70

Natural (unregulated) streamflow

EXPLANATION

Simulated streamflow

Figure 17. Simulated and natural (unregulated) monthly streamflow for Long Creek at the western crossing of the 

international boundary (05NA003).

A. Monthly mean streamflow 

Monthly mean streamflow,

in cubic meters per second

0

2

4

6

8

Month

Monthly standard deviation of streamflow,

in cubic meters per second

Jan. Feb. Mar. Apr. May June July Aug. Sept. Oct. Nov. Dec.

0

2

4

6

8

10

12

B. Monthly standard deviation of streamflow

Natural (unregulated) streamflow

EXPLANATION

Simulated streamflow

Natural (unregulated) streamflow

EXPLANATION

Simulated streamflow

Figure 18. Simulated and natural (unregulated) monthly streamflow for Long Creek at the western crossing of the international 

boundary (05NA003). A, monthly mean streamflow. B, monthly standard deviation of streamflow.
32 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

Table 9. Lags applied to monthly streamflow routing model.

[Sask., Saskatchewan; N. Dak., North Dakota]

Station/subbasin 

identifier Station name State/province

Fraction of current 

monthly streamflow 

Fraction of previous 

monthly streamflow

05ND004 Moose Mountain Creek near Oxbow Sask. 1.0 0

05NA003 Long Creek at western crossing of international boundary Sask. 1.0 0

05113600 Long Creek near Noonan N. Dak. 1.0 0

05NB001 Long Creek near Estevan Sask. 0.8 0.2

05NB036 Souris River below Rafferty Reservoir Sask. 0.8 0.2

05114000 Souris River near Sherwood N. Dak. 0.95 0.05

05116000 Souris River near Foxholm N. Dak. 1.0 0

05116500 Des Lacs River at Foxholm N. Dak. 1.0 0

05117500 Souris River above Minot N. Dak. 0.9 0.1

05120000 Souris River near Verendrye N. Dak. 0.8 0.2

05120500 Wintering River near Karlsruhe N. Dak. 0.8 0.2

05122000 Souris River near Bantry N. Dak. 0.75 0.25

05123400 Willow Creek near Willow City N. Dak. 0.75 0.25

05123510 Deep River near Upham N. Dak. 0.75 0.25

05124000 Souris River near Westhope N. Dak. 1.0 0

words, historic ratios during a specific year and month were 

applied to the same year and month for the simulated values. 

Once all the monthly subbasin flows had been subdivided into 

three-per-month flows, they were routed through the Souris 

River Basin in a similar manner as the monthly flows, with 

additional consideration taken for streamflow loss to evaporation and aquifer recharge (table 10). The three-per-month flows 

were calibrated such that simulated 10-day means and standard 

deviations compared well with the natural flow data provided 

by the USACE (U.S. Army Corps of Engineers, 2013).

At the three-per-month time step, it was easier to catch 

small losses in streamflow assumed to be contributed to evaporation and (or) aquifer recharge. Subbasins with lagged flows 

greater than 1 day had noticeable streamflow losses, which 

increased as the lag increased (table 10). This may be the result 

of increased exposure of streamflow to surface evaporation 

as well as the result of streamflow loss to nearby/underlying 

aquifers. Specifically, Pusc (1994) noted that the Souris River 

can lose anywhere from 0.13 to 0.52 m3

/s to the Minot aquifer 

(fig. 14), with a potential for losing greater amounts (Pettyjohn, 

1967; Bradley, 1963). The Minot aquifer is within the Souris 

River near Verendrye, N. Dak. (05120000) subbasin (fig. 14), 

but losses to the Minot aquifer were too small for the three-permonth routing model to catch; however, Pusc (1994) focused 

mainly on water losses to the Minot aquifer and there are 

several aquifers that lie downstream, which may contribute to 

streamflow loss. These aquifers include, but are not limited to, 

Sundre, Voltaire, Karlsruhe, and New Rockford (fig. 14). 

In addition to considering the effects of evaporation and 

aquifer recharge in the three-per-month streamflow routing 

model, the effects of the J. Clark Salyer NWR also resulted in 

modifications to the routing model. Inclusion of the J. Clark 

Salyer NWR resulted in a reduction of peak flows during the 

spring and increased flows during the summer at the Souris 

River near Westhope, N. Dak., (05124000) gaging station 

because of water being retained in the spring for management 

of vegetation and then slowly released throughout the summer 

months (U.S. Army Corps of Engineers, 2013). Although the 

amount of water retained for management of vegetation in the 

spring and the amount released throughout the summer varied 

each year depending on management needs and water availability, an average amount retained and released was determined through trial and error such that simulated 10-day means 

and standard deviations compared well with the natural 10-day 

historical flows provided by the USACE (U.S. Army Corps 

of Engineers, 2013). To account for the effects of the J. Clark 

Salyer NWR in the three-per-month streamflow routing model, 

peak flows were reduced in the spring by retaining 15 percent 

of the natural streamflow during March, April, and May. Then, 

during the summer months (June, July, August, September), 

some of the retained flow was released such that 8 percent of 

the average 10-day flow between March and May was released 

for each 10-day period between June and September, resulting in slightly higher flows at Westhope during the summer 

months. A list of basin-specific lags and streamflow loss coefficients is provided in table 10. 

Routed 10-day simulated flows compared well with 

natural flows provided by the USACE, with 10-day means and 

standard deviations showing little bias and following similar 

trends in 10-day variation (figs. 21 and 22); hence, both routed 

10-day and monthly natural flows matched the statistical properties of historical natural flow quite well.
Water-Balance Model for Estimating Natural Streamflow 33

0

10

20

30

40

10

20

30

40

50

60

70

Monthly mean streamflow,

in cubic meters per second

Month

Monthly mean streamflow,

in cubic meters per second

Jan. Feb. Mar. Apr. May June July Aug. Sept. Oct. Nov. Dec.

0

A. Souris River near Sherwood, North Dakota (05114000)

B. Souris River near Westhope, North Dakota (05124000)

Natural (unregulated) streamflow

EXPLANATION

Simulated streamflow

Natural (unregulated) streamflow

EXPLANATION

Simulated streamflow

Figure 19. Means of simulated and natural (unregulated) monthly streamflow. A, the Souris River near 

Sherwood, North Dakota (05114000) streamflow-gaging station. B, the Souris River near Westhope, North Dakota 

(05124000) streamflow-gaging station.
34 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

10

0

20

30

40

50

10

20

30

40

50

60

70

80

Monthly standard deviation of streamflow,

in cubic meters per second

Month

Monthly standard deviation of streamflow,

in cubic meters per second

Jan. Feb. Mar. Apr. May June July Aug. Sept. Oct. Nov. Dec.

0

A. Souris River near Sherwood, North Dakota (05114000)

B. Souris River near Westhope, North Dakota (05124000)

Natural (unregulated) streamflow

EXPLANATION

Simulated streamflow

Natural (unregulated) streamflow

EXPLANATION

Simulated streamflow

Figure 20. Standard deviations of simulated and natural (unregulated) monthly streamflow. A, the Souris River 

near Sherwood, North Dakota (05114000) streamflow-gaging station. B, Souris River near Westhope, North Dakota 

(05124000) streamflow-gaging station.
Water-Balance Model for Estimating Natural Streamflow 35

Table 10. Lags applied to 10-day streamflow routing model.

[Sask., Saskatchewan; N. Dak., North Dakota; NA, not available due to lack of reconstructed streamflow data]

Station/subbasin 

identifier Station name State/province

Fraction of 

current 10-day 

streamflow 

Fraction of 

previous 10-day 

streamflow

Streamflow loss 

(percent)

05ND004 Moose Mountain Creek near Oxbow Sask. 1.0 0 0

05NA003 Long Creek at western crossing of international 

boundary

Sask. 1.0 0 0

05113600 Long Creek near Noonan N. Dak. 0.80 0.20 5

05NB001 Long Creek near Estevan Sask. NA NA NA

05NB036 Souris River below Rafferty Reservoir Sask. 0.80 0.20 5

05114000 Souris River near Sherwood N. Dak. 0.90 0.10 0

05116000 Souris River near Foxholm N. Dak. 1.0 0 0

05116500 Des Lacs River at Foxholm N. Dak. 1.0 0 0

05117500 Souris River above Minot N. Dak. 0.90 0.10 0

05120000 Souris River near Verendrye N. Dak. 0.65 0.35 7

05120500 Wintering River near Karlsruhe N. Dak. 0.65 0.35 7

05122000 Souris River near Bantry N. Dak. 0.60 0.40 (1) 

05123400 Willow Creek near Willow City N. Dak. 0.60 0.40 (1) 

05123510 Deep River near Upham N. Dak. 0.70 0.30 (1) 

05124000 Souris River near Westhope N. Dak. 1.0 0 0

1

Streamflow loss is dependent on the amount of flow retained in the J. Clark Salyer National Wildlife Refuge versus the amount of flow released during the 

summer months.
36 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

0

20

40

60

80

0

20

40

60

80

10-day mean streamflow, in cubic meters per second 10-day mean streamflow, in cubic meters per second

Month

Jan. Feb. Mar. Apr. May June July Aug. Sept. Oct. Nov. Dec.

A. Souris River near Sherwood, North Dakota (05114000)

B. Souris River near Westhope, North Dakota (05124000)

Natural (unregulated) streamflow

EXPLANATION

Simulated streamflow

Natural (unregulated) streamflow

EXPLANATION

Simulated streamflow

Figure 21. Means of simulated and natural (unregulated) 10-day streamflow. A, the Souris River near Sherwood, North 

Dakota (05114000) streamflow-gaging station. B, Souris River near Westhope, North Dakota (05124000) streamflow-gaging 

station.
Water-Balance Model for Estimating Natural Streamflow 37

0

20

40

60

80

100

120

0

20

40

60

80

100

120

10-day standard deviation of streamflow,

in cubic meters per second

Month

10-day standard deviation of streamflow,

in cubic meters per second

Jan. Feb. Mar. Apr. May June July Aug. Sept. Oct. Nov. Dec.

A. Souris River near Sherwood, North Dakota (05114000)

B. Souris River near Westhope, North Dakota (05124000)

Natural (unregulated) streamflow

EXPLANATION

Simulated streamflow

Natural (unregulated) streamflow

EXPLANATION

Simulated streamflow

Figure 22. Standard deviations of simulated and natural (unregulated) 10-day streamflow. A, the Souris River near 

Sherwood, North Dakota (05114000) streamflow-gaging station. B, the Souris River near Westhope, North Dakota (05124000) 

streamflow-gaging station.
38 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

Stochastic Natural Streamflow Model

When evaluating future flood risk for the Souris River 

Basin, revising reservoir operating rules, or designing levees 

or other hydraulic structures, it is critical to have a model that 

can simulate future flows that realistically mimic natural flows 

and accurately represent the probabilities of future extreme 

events. The stochastic natural streamflow model is intended to 

simulate future conditions well beyond the envelope of those 

seen in the historical record and assign probabilities that can 

be used to evaluate risk and complete economic cost/benefit 

analyses. In this section, the stochastic climate simulation 

and water-balance models, described in previous sections, are 

combined to simulate future natural flows, which in turn are 

used to compute flow-frequency curves (frequency distributions). Such distributions are used extensively for designing 

hydraulic structures and completing cost/benefit analyses; 

however, in the standard application of these methods it is 

assumed that both climate and streamflow are stationary (that 

is, the frequency distribution remains constant through time). 

Long historical streamflow records are used to compute an 

empirical frequency distribution, which in turn is modeled 

using a theoretical statistical distribution such as the log-Pearson type III distribution (Interagency Advisory Committee on 

Water Data, 1982); however, for cases such as the Souris River 

Basin where (1) climate is not stationary and (2) extreme 

hydroclimatic variability and persistence make determining 

the theoretical frequency distribution difficult even with a long 

historical record, it is necessary to use a different approach 

such as the one described in this section.

Methods

To combine the stochastic climate simulation model 

(SCSM) and the WBM, the group-averaged seasonal data from 

the SCSM (eq. 1) needed to be converted into monthly data for 

specific locations that could be used by the WBM. To ensure 

that the simulated monthly data realistically represented the 

statistical properties (means, variances, serial correlations, and 

cross-correlations) of actual data, a simple conditional resampling approach was used: 

PET y R y PET y P y

R y P y L

L m L m

PET

G S L m

L m

P

G S

, , , ,

, ,

( ) ( ) ( ); ( )

( ) ( ); ,

=

= = 1   , ; 4 1 m = , ,4 (10)

where

PET y L m, ( ) is the potential evapotranspiration for the 

meteorological station L, month m, and 

year y;

R y L m

PET

, ( ) is the ratio computed from historical data 

for potential evapotranspiration PET, 

meteorological station L, month m, and 

year y;

P y L m, ( ) is the precipitation for the meteorological 

station L, month m, and year y; and

R y L m

P

, ( ) is the ratio computed from historical data for 

precipitation P, meteorological station L, 

month m, and year y.

where L designates the particular meteorological station 

locations in group G and m designates the particular months 

in season S. For each year, there were 192 ratios (16 for each 

of the 12 combinations of station group and season) for PET 

and 192 ratios for precipitation. For each year of the historical record from y=1946 to y=2011, the ratios were computed. 

The starting year 1946 was used because ratios needed to be 

matched with the ratios used for disaggregating simulated 

flows from a monthly to 10-day time step, and the latter ratios 

were not available before 1946. For a given simulation year, 

the group-averaged seasonal values for PET and precipitation 

generated using the SCSM were multiplied by the ratios for a 

randomly sampled (with replacement) historical year to obtain 

simulated values of PET and precipitation for each month and 

station location:

PET R y PET P

R y P L

L m L m

PET

G S L m

L m

P

G S

, , , ,

, ,

(Y) ( *) (Y); (Y)

( *) (Y);

=

= = 1 4 , ,   ; , m = 1 4, (11)

where Y designates a simulation year and y* a randomly 

sampled historical year. In accordance with the dry (1912–69) 

and wet (1970–2011) climate states described previously for 

the SCSM, for the dry climate state, y* is randomly sampled 

from historical years 1946–69, and for the wet climate state, 

y* is randomly sampled from historical years 1970–2011. 

Simulated monthly temperature was back-calculated from 

simulated monthly PET using separate equations for each 

combination of month and station location. Simulated values 

for monthly PET, precipitation, and temperature for each 

meteorological station location and month were then interpolated to the 8x8-km grid cells using the same loess procedure 

described previously for the WBM calibration and used to 

obtain simulated monthly streamflow values for each subbasin. 

The monthly streamflow values for each subbasin were disaggregated to three values per month using the same randomly 

sampled historical year y* that was used to disaggregate the 

climatic inputs, and the 10-day streamflow routing model 

was used to simulate streamflow for each gaging station. 

The reason for using the same randomly sampled year is that 

the relative distribution of monthly flow among the 10-day 

subintervals for a particular year should be determined by 

climatic conditions (timing of snowmelt, for example) for that 

year. Although the climate simulations consisted of “climate 

years” November–October (November and December from 

the previous calendar year combined with January–October of 

the current calendar year), the WBM simulations were done 

using calendar years (January–December). Because November 

and December are winter months with minimal runoff, the difference between using calendar year y*or y*–1 ratios for those 

months to disaggregate monthly flows was expected to be negligible; furthermore, because in either case (climate years or 

calendar years) the break between years happens in late fall or 
Stochastic Natural Streamflow Model 39

winter, discontinuity from the end of one year to the beginning 

of the next year caused by randomly sampling 2 separate years 

to disaggregate climate variables or streamflow was minimal.

Stochastic Simulation Results

As discussed in the “Stochastic Climate Model for Simulation of Precipitation, Temperature, and Potential Evapotranspiration” section of this report, the Souris River Basin can be 

characterized by wet or dry equilibrium climate states. Soilmoisture storage can potentially cause a long lag-time between 

the onset of a new climate state and the eventual onset of the 

new streamflow equilibrium. To evaluate flood risk in this light, 

100 independent realizations, or traces, of streamflow data 

were generated, each starting in a wet equilibrium state that 

continues for 50 years before transitioning to a dry state that 

continues for an additional 50 years. It is not necessary at the 

outset to specify exactly when this transition may take place in 

the future or even to specify a starting year for the simulations. 

To avoid dependence of the WBM simulations on initial conditions (such as starting soil-moisture storage), an initial 10-year 

“burn-in” period was discarded at the beginning of each trace. 

Example simulated traces are illustrated for three of 

the gaging stations (table 6)—Souris River below Rafferty 

Reservoir (05NB036), Souris River above Minot, N. Dak. 

(05117500), and Wintering River near Karlsruhe, N. Dak. 

(05120500) (figs. 1 and 14; hereafter referred to as the “Rafferty,” “Minot,” and “Wintering” gaging stations, respectively). 

For comparison with the actual data, historical natural streamflows for 1946–2011 are shown in figure 23. The red lines in 

these graphs show the median, 90th percentile, and maximum 

recorded values for the annual maximum 10-day natural 

streamflow during 1946–2011. These graphs highlight the 

extreme interannual streamflow variability of the Souris River 

Basin; for example, for the Minot gaging station, the annual 

maxima are near zero for many years and for a “normal” year 

(represented by the median), the maximum is about 30 m3/s. 

To put this value in perspective, 30 m3

/s of flow for a 10-day 

period for this site corresponds to only about 1 mm of runoff 

for the upstream drainage area. At the upper extreme, in 2011 

the maximum 10-day flow was about 500 m3/s, which corresponds to about 16 mm of runoff. During most years, runoff for 

the Souris River Basin is only a small fraction of the precipitation that falls in the basin and is much less than the soil-moisture storage capacity; thus, streamflow is extremely variable 

and persistent (characterized by a tendency for clusters of 

low-flow and high-flow years). The same properties are evident 

in the historical flows for the other two gaging stations as well. 

One of the 100 simulated streamflow traces for the three 

selected gaging stations is shown in figure 24. This trace is 

similar overall to the natural streamflow record (fig. 23) (but 

with the wet period at the beginning rather than the end). For 

the Rafferty gaging station, during the wet simulation period 

there was one simulation year that exceeded the 2011 recorded 

maximum, and 10 percent (5 out of 50) of the annual maxima 

exceeded the 90th percentile of the recorded maxima. For the 

other two gaging stations and the wet simulation period, none 

of the simulated annual maxima exceeded the value from 2011 

but several simulated maxima exceeded the 90th percentile 

line. There was a stark contrast for this trace between the simulated flows for the wet and dry periods. For the dry period, only 

one of the simulated maxima for both the Rafferty and Minot 

gaging stations and three for the Wintering gaging station 

exceeded the 90th percentile line.

Another simulated trace is shown in figure 25. This 

trace differs from the previous example in that there are more 

extreme simulated high flows during the wet period. For the 

Rafferty and Minot gaging stations, there are two simulation 

years during the wet period for which the annual maxima 

exceeded the value from 2011. In addition, for the Rafferty and 

Minot gages the simulated maximum for the second year of the 

dry period exceeded the 2011 value, which indicates that high 

soil moisture conditions at the end of the wet period can persist 

for several years after the transition to drier conditions; however, as shown by the high simulated value for Rafferty near 

the end of the dry period, occasional extremes similar to 2011 

can take place even after the lower streamflow equilibrium is 

reached. 

The set of 100 independent streamflow traces were used 

to estimate equilibrium frequency distributions for the wet 

and dry climate states. The first 50 years of each trace (5,000 

simulation years for the 100 traces) are assumed to represent 

the equilibrium distribution for the wet climate state, and the 

last 40 years of each trace (4,000 simulation years) are assumed 

to represent the equilibrium distribution for the dry climate 

state. Years 51–60 are assumed to represent a transition period 

from the wet to the dry state before equilibrium was reached 

for the dry state. The probability of exceeding a given value for 

a specified streamflow variable (such as the annual maximum 

10-day flow or the annual total flow volume) was estimated 

using the percent of the simulation years for which the value 

was exceeded.

The equilibrium frequency distributions for annual maximum 10-day mean streamflow and annual total streamflow volume for the Rafferty gaging station are shown in figure 26. To 

allow comparison with the historical data, empirical frequency 

distributions for the historical natural flow data for 1946–69 

(representing the dry climate state) and 1970–2011 (representing the wet climate state) are plotted along with the curves 

computed using the simulated data. The Rafferty gaging station 

location is important because it represents natural inflows to 

Rafferty Reservoir, which is the major flood control structure 

for downstream communities. The capacity of Rafferty Dam 

at its maximum storage capacity is 633 million cubic meters 

(International Joint Commission, 1989), which is shown by the 

horizontal line in figure 26B. In 2011, total inflow to Rafferty 

Reservoir was almost 1.2 billion cubic meters, twice its capacity. Much of the inflow came during a span of a few weeks in 

May and June of that year, overwhelming the reservoir and 

causing a large uncontrolled spill. The equilibrium frequency 

curve for annual total volume for the wet climate state (the 
40 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

Souris River above Minot, North Dakota (05117500)

0

100

200

300

400

500

600

Souris River below Rafferty Reservoir, Saskatchewan (05NB036)

0

50

100

150

200

250

Streamflow, in cubic meters per second

Wintering River near Karlsruhe, North Dakota (05120500)

0

5

10

15

20

25

30

35

40

1945 1950 1955 1960 1965 1970 1975 1980

Year

1985 1990 1995 2000 2005 2010 2015

Q100

Q90

Q50

Q100

Q90

Q50

Q100

Q90

Q50

Figure 23. Natural (unregulated) 10-day mean streamflow for 1946–2011 for the Souris River below Rafferty Reservoir, 

Saskatchewan (05NB036), Souris River above Minot, North Dakota (05117500), and Wintering River near Karlsruhe, 

North Dakota (05120500), streamflow-gaging stations as compared to the median (Q50), 90th percentile (Q90), and 

maximum (Q100) recorded values of annual maximum 10-day streamflow.
Stochastic Natural Streamflow Model 41

Souris River above Minot, North Dakota (05117500)

0

100

200

300

400

500

600

Souris River below Rafferty Reservoir, Saskatchewan (05NB036)

0

50

100

150

200

250

Streamflow, in cubic meters per second

Wintering River near Karlsruhe, North Dakota (05120500)

0

5

10

15

20

25

30

35

40

0 10 20 30 40 50

Simulation year

60 70 80 90 100

Q100

Wet period Dry period

Wet period Dry period

Wet period Dry period

Q90

Q50

Q100

Q90

Q50

Q100

Q90

Q50

Figure 24. Simulated 10-day mean streamflow for 100-year simulation period (50 years of wet climate state followed 

by 50 years of dry climate state) for the Souris River below Rafferty Reservoir, Saskatchewan (05NB036), Souris River 

above Minot, North Dakota (05117500), and Wintering River near Karlsruhe, North Dakota (05120500), streamflow-gaging 

stations (trace 10 out of 100) as compared to the median (Q50), 90th percentile (Q90), and maximum (Q100) recorded 

values of annual maximum 10-day streamflow.
42 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

Souris River above Minot, North Dakota (05117500)

0

100

200

300

400

500

600

Souris River below Rafferty Reservoir, Saskatchewan (05NB036)

0

50

100

150

200

250

Streamflow, in cubic meters per second

Wintering River near Karlsruhe, North Dakota (05120500)

0

5

10

15

20

25

30

35

40

0 10 20 30 40 50 60 70 80 90 100

Q100

Q90

Q50

Q100

Q90

Q50

Q100

Q90

Q50

Simulation year

Wet period Dry period

Wet period Dry period

Wet period Dry period

Figure 25. Simulated 10-day mean streamflow for 100-year simulation period (50 years of wet climate state followed 

by 50 years of dry climate state) for the Souris River below Rafferty Reservoir, Saskatchewan (05NB036), Souris River 

above Minot, North Dakota (05117500), and Wintering River near Karlsruhe, North Dakota (05120500), streamflow-gaging 

stations (trace 4 out of 100) as compared to the median (Q50), 90th percentile (Q90), and maximum (Q100) recorded 

values of annual maximum 10-day streamflow.
Stochastic Natural Streamflow Model 43

0.1

1

10

100

1,000

0.2

2

20

200

0.5

5

50

500

Streamflow, in cubic meters per second

A. Annual maximum 10-day mean streamflow

0.1

1

10

100

1,000

0.2

2

20

200

2,000

0.5

5

50

500

5,000

Volume, in million cubic meters

99.8 99.5 99 95 90 80 70 60 50 40 30 20 10 5 1 0.5 0.2

Exceedance probability, in percent

B. Annual total streamflow volume

Wet (1970−2011)

Dry (1912−69)

Natural

(unregulated)

Modeled

EXPLANATION

Natural

(unregulated)

Modeled

EXPLANATION

Volume of Rafferty Reservoir at maximum storage capacity

Dry (1912−69)

Wet (1970−2011)

Figure 26. Equilibrium frequency distributions for the Souris River below Rafferty Reservoir, Saskatchewan (05NB036), 

streamflow-gaging station for wet (1970–2011) and dry (1912–1969) climate states. A, annual maximum 10-day streamflow; 

and B, annual total streamflow volume. 
44 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

blue curve in figure 26B) indicates only about a 0.2 percent per 

year (1 in 500) chance of exceeding the 2011 inflow volume; 

however, it indicates a much more likely chance—about 3 

percent per year—of exceeding the 633 million cubic meter 

capacity. In contrast, the dry equilibrium curve indicates only 

about a 0.5-percent chance of exceeding the 633 million cubic 

meter capacity; thus, it is about six times more likely to exceed 

633 million cubic meters of inflow for the wet state compared 

to the dry state. Considering that the wet state could continue 

for several more decades, this high flood risk could have major 

implications for downstream flood control.

The equilibrium frequency distributions for the Minot 

gaging station are shown in figure 27. The maximum 10-day 

natural flow for this site in 2011 was about 500 m3/s, which 

falls just below the equilibrium frequency curve for the wet 

state, indicating about a 3-percent per year chance of exceedance. In comparison, the dry equilibrium curve indicates less 

than a 0.5-percent per year chance of exceeding the same 

value. The annual total natural flow volume for the Minot 

gaging station was about 2.4 billion cubic meters in 2011, 

which, similar to the Rafferty gaging station, is estimated to 

have about a 0.5-percent per year chance of exceedance for the 

wet state. For the dry state, there is an extremely small chance 

(less than 0.1 percent) of exceeding the 2011 flow volume.

The equilibrium frequency distributions for the Westhope 

gaging station are shown in figure 28. This gage is at the 

international border between North Dakota and Manitoba. 

The results for this gage are similar to the previous examples; 

thus, it seems that the wet and dry climate states have a similar 

effect on streamflow for the entire main stem Souris River 

upstream from Westhope, which reinforces the conclusion that 

the shifting climate states represent a large regional phenomenon encompassing the Souris River Basin and other areas of 

Saskatchewan, North Dakota, and Manitoba.

1

10

100

1,000

2

20

200

5

50

500

1

10

100

1,000

2

20

200

2,000

5

50

500

5,000

Streamflow, in cubic meters per second Volume, in million cubic meters

99.8 99.5 99 95 90 80 70 60 50 40 30 20 10 5 1 0.5 0.2

Exceedance probability, in percent

Wet (1970−2011)

Dry (1912−69)

Dry (1912−69)

Wet (1970−2011)

A. Annual maximum 10-day mean streamflow

B. Annual total streamflow volume

Natural

(unregulated)

Modeled

EXPLANATION

Natural

(unregulated)

Modeled

EXPLANATION

Figure 27. Equilibrium frequency distributions for the Souris River above Minot, North Dakota (05117500), streamflow-gaging 

station for wet (1970–2011) and dry (1912–69) climate states. A, annual maximum 10-day mean streamflow; and B, annual total 

streamflow volume.
Summary 45

Streamflow, in cubic meters per second Volume, in million cubic meters

99.8 99.5 99 95 90 80 70 60 50 40 30 20 10 5 1 0.5 0.2

Exceedance probability, in percent

Wet (1970−2011)

Dry (1912−69)

Dry (1912−69)

Wet (1970−2011)

A. Annual maximum 10-day mean streamflow

B. Annual total streamflow volume

Natural

(unregulated)

Modeled

EXPLANATION

Natural

(unregulated)

Modeled

EXPLANATION

1

10

100

1,000

2

20

200

2,000

5

50

500

5,000

1

10

100

1,000

2

20

200

2,000

5

50

500

5,000

Figure 28. Equilibrium frequency distributions for the Souris River near Westhope, North Dakota (05124000), streamflowgaging station for wet (1970–2011) and dry (1912–69) climate states. A, annual maximum 10-day mean streamflow; and B, 

annual total streamflow volume.

Summary

The Souris River Basin is a 61,000 square-kilometer basin 

in the Provinces of Saskatchewan and Manitoba and the State 

of North Dakota. In May and June of 2011, record-setting 

rains were seen in the headwater areas of the basin. Emergency 

spillways of major reservoirs were discharging at or near full 

capacity resulting in extensive flooding to numerous downstream communities. As a result, a Souris River Basin Task 

Force was appointed to develop a plan of study for mitigating 

future damages from extreme events such as the 2011 flood. 

One of the recommendations of the task force was to develop 

a model for simulating future streamflow that could determine 

the probabilities of future floods or droughts. In response to this 

recommendation, the U.S. Geological Survey, in cooperation 

with the North Dakota State Water Commission, developed a 

stochastic model for simulating Souris River Basin precipitation, evapotranspiration, and natural (unregulated) streamflow. 

Simulations from the model can be used in future studies to 

simulate regulated streamflow, design levees, and other structures, and complete economic cost/benefit analyses. 

Long-term climatic variability was analyzed using treering chronologies and meteorological data from United States 

and Canadian weather stations in a large region of southcentral Canada and the north-central United States centered 

on the Souris River Basin. Five statistically distinct clusters of 

weather stations were identified based on seasonal precipitation 

patterns, resulting in five distinct climate regions, four of which 

overlapped the Souris River Basin. Using multiresolution 

decomposition of tree-ring chronologies along with multiple 

regression equations to predict seasonal precipitation based on 

tree-ring wavelets with wavelengths of 16, 32, and 64 years, it 
46 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

was determined that most of the multidecadal variability of seasonal precipitation for the climate regions could be explained 

by the tree-ring chronologies. The tree-ring chronologies were 

used to hindcast precipitation to the early 1700s and compare 

recent (since the late 1800s) wet and dry conditions to earlier extreme conditions. It was shown that several extreme 

wet periods (as extreme as or more extreme than the wettest 

recorded period) and several extreme dry periods (as extreme 

as or more extreme than the driest recorded period) took place 

prior to the historical record. The extreme wet and dry periods 

were verified based on historical information and comparison 

with other paleoclimatic studies.

A stochastic climate model for simulation of precipitation, 

temperature, and potential evapotranspiration for the Souris 

River Basin was developed using the recorded meteorological data and the tree-ring-extended precipitation records. The 

model was developed first for seasonal data (three 4-month 

seasons: November–February, March–June, and July–October), 

and then the seasonal data were disaggregated to a monthly 

time step. A significant climate transition was seen around 

1970, with 1912–69 representing a dry climate state and 1970–

2011 representing a wet climate state. Although there were 

some distinct subpatterns within the basin, the predominant 

differences between the two states were higher spring through 

early fall precipitation and higher spring potential evapotranspiration for the wet compared to the dry state. Verification of 

the stochastic climate simulation model showed that it was able 

to reproduce both the short-term (annual to multidecadal) variability of the historical data and the long-term (multidecadal to 

century scale) variability of the tree-ring-extended data.

A water-balance model was developed for simulating 

monthly natural (unregulated) streamflow based on precipitation, temperature, and potential evapotranspiration. The model 

was calibrated using streamflow data from the U.S. Geological 

Survey and Environment Canada, as well as natural (unregulated) streamflow data from the U.S. Army Corps of Engineers. 

Streamflow was simulated for 15 streamflow-gaging stations. 

The water-balance model was designed to estimate monthly 

runoff from each of 1,103 8x8-kilometer grid cells covering the 

Souris River Basin based on moisture inputs (precipitation), 

moisture outputs (evapotranspiration and runoff), and soilmoisture storage. Static inputs include soil-moisture storage 

capacity and soil permeability; and dynamic inputs include 

monthly precipitation, temperature, and potential evapotranspiration interpolated to each grid cell from weather station 

data using a locally weighted smoothing algorithm. Monthly 

runoff for the grid cells was aggregated to obtain runoff for 

each of 15 subbasins, and a simple streamflow routing model 

was used to simulate monthly runoff for streamflow-gaging 

stations with more than one upstream subbasin. After calibrating the model for a monthly time step, monthly runoff for each 

subbasin was disaggregated into three values per month, or 

an approximately 10-day time step, and a separate streamflow 

routing model was developed for simulating 10-day streamflow 

for downstream streamflow-gaging stations. Water-balance 

model results indicated good agreement between simulated 

and natural (unregulated) streamflow for both the monthly and 

10-day time steps. Correlation coefficients between simulated 

and natural (unregulated) flows were generally high and the 

seasonal means and standard deviations of the simulated flows 

closely matched the means and standard deviations of the natural (unregulated) flows.

The stochastic climate simulation model for precipitation, temperature, and potential evapotranspiration was 

combined with the water-balance model to simulate potential 

future sequences of 10-day mean streamflow for each of the 

streamflow-gaging station locations. Flood risk, as determined 

by equilibrium flow-frequency distributions for both the dry 

(1912–69) and wet (1970–2011) climate states, is considerably 

higher for the wet state compared to the dry state; for example, 

for the wet climate state, annual inflow volume for Rafferty 

Reservoir (the major flood control structure for downstream 

communities) has about a 3-percent chance per year of exceeding 633 million cubic meters, the maximum storage capacity of 

the reservoir. In contrast, for the dry climate state, the annual 

inflow volume has only about a 0.5-percent chance per year 

of exceeding 633 million cubic meters. Future flood risk will 

remain high until the wet climate state ends, and for several 

years after that, because there may be a long lag-time between 

the return of drier conditions and the onset of a lower soilmoisture storage equilibrium. 

References Cited

Agriculture and Agri-Food Canada, 2013, Soil Landscapes 

of Canada, version 3.2: Government of Canada, accessed 

October 21, 2014, at http://sis.agr.gc.ca/cansis/nsdb/slc/

v3.2/index.html.

Agriculture and Agri-Food Canada, 2014, Detailed Soil Survey (DSS) compilations: Government of Canada, accessed 

October 21, 2014, at http://sis.agr.gc.ca/cansis/nsdb/dss/v3/

index.html.

Ault, T.R., and St. George, Scott., 2010, The magnitude of 

decadal and multidecadal variability in North American 

precipitation: Journal of Climate, v. 23, p. 842–850.

Bluemle, J.P., 1996, From the state geologist: North Dakota 

Geological Survey Newsletter, v. 23, no. 1, p. 1–2.

Bradley, Edward., 1963, Relation of surface and ground water 

in the Souris River Valley near Minot, North Dakota: U.S. 

Geological Survey Open-File Report 62–10, 12 p.

Bunn, A.G., Korpela, Mikko; Biondi, Franco; Campelo, Filipe; 

Mérian, Pierre; Mudelsee, Manfred; Qeadan, Fares; Schulz, 

Michael; and Zang, Christian, 2014, dplR—Dendrochronology program library in R: R Foundation for Statistical 

Computing, R package version 1.6.0.
References Cited 47

Cleveland, W.S., Grosse, Eric., and Shyu, W.M., 1992, Local 

regression model, chap. 8 of Chambers, J.M., and Hastie, 

T.J., eds., Statistical models in S: Pacific Grove, California, 

Wadsworth & Brooks/Cole, p., 309–376.

Environment Canada, 2014, Adjusted and Homogenized Canadian Climate Data (AHCCD): Environment and Climate 

Change Canada, accessed July 15, 2015, at https://ec.gc.ca/

dccha-ahccd/default.asp?lang=En&n=B1F8423A-1.

Esri, 2014 ArcMap 10.2.2 Hydro Tools: Esri, accessed October 9, 2014, at http://blogs.esri.com/esri/arcgis/2011/10/12/.

Garbrecht, J.D., and Rossel, F.E., 2002, Decade-scale precipitation increase in Great Plains at end of 20th century: 

Journal of Hydrologic Engineering, v. 7, no. 1, p. 64–75.

Government of Canada, 2015a, PFRA sub-basins of the 

AAFC Watersheds Project—2013: Government of Canada, 

accessed June 11, 2015, at http://open.canada.ca/data/en/

dataset/4f3c7d6d-e018-4a69-a6cf-a4c327572b24.

Government of Canada, 2014, Water office: Government of 

Canada, accessed November 19, 2014, at http://wateroffice.

ec.gc.ca/.

Government of Canada, 2015b, Total gross drainage areas 

of the AAFC Watersheds Project—2013: Government of 

Canada, accessed June 11, 2015, at http://open.canada.ca/

data/en/dataset/67c8352d-d362-43dc-9255-21e2b0cf466c.

Gray, S.T., and McCabe, G.J., 2010, A combined water balance and tree ring approach to understanding the potential 

hydrologic effects of climate change in the central Rocky 

Mountain region: Water Resources Research, v. 46, no. 5, 

W05513, accessed August 20, 2015, at http://onlinelibrary.

wiley.com/doi/10.1029/2008WR007650/full.

Hamon, W.R., 1961, Estimating potential evapotranspiration: 

Proceedings of the American Society of Civil Engineers, 

Journal of the Hydraulic Division, v. 87, no. 3, p. 107–120.

Helsel, D.R., and Hirsch, R.M., 2002, Statistical methods in 

water resources: U.S. Geological Survey Techniques of 

Water-Resources Investigations, book 4, chap. A3, 522 p.

Hirsch, R.M., and Ryberg, K.R., 2012, Has the magnitude of 

floods across the USA changed with global CO2 levels?: 

Hydrological Science Journal, v. 57, no. 1, p. 1–9.

Hoerling, Martin; Eischeid, Jon; Esterling, David; Peterson, Thomas; and Webb, Robert, 2010, Understanding 

and explaining hydro-climatic variations of Devils Lake: 

National Oceanic and Atmospheric Administration, Climate 

Assessment Report, 20 p.

Huffman, R.L., Fangmeier, D.D., Elliot, W.J., Workman, 

S.R., and Schwab, G.O., 2011, Soil and water conservation 

engineering (6th ed.): St. Joseph, Mich., American Society 

of Agricultural Engineers, 523 p.

Interagency Advisory Committee on Water Data, 1982, Guidelines for determining flood flow frequency: Reston, Va., 

U.S. Geological Survey Hydrology Subcommittee Bulletin 

17B, 28 p. plus appendixes. [Also available at http://water.

usgs.gov/osw/bulletin17b/bulletin_17B.html.]

International Joint Commission, 1989, Boundary waters—

Agreement between the Government of Canada and the 

Government of the United States of America for water 

supply and flood control in the Souris River Basin (with 

Annexes and Canada/Sasketchewan agreement): Ottawa, 

Canada, Queen’s Printer for Canada, 36 p.

International Souris River Board, 2013, Plan of study—For 

the review of the operating plan contained in Annex A of 

the 1989 international agreement between the Government 

of Canada and the Government of the United States of 

America: Souris River Basin Task Force, 197 p.

Johnson, R.C., Mercier, T.J., and Brownfield, M.E., 2014, 

Spatial and stratigraphic distribution of water in oil shale of 

the Green River Formation using Fischer Assay, Piceance 

Basin, northwestern Colorado: U.S. Geological Survey 

Open-File Report 2014–1059, 57 p., accessed April 14, 

2014, at http://dx.doi.org/10.3133/ofr20141059.

Lu, Jianbiao; Sun, Ge; McNulty, S.G., and Amatya, D.M., 

2005, A comparison of six potential evapotranspiration 

methods for regional use in the southeastern United States: 

Journal of the American Water Resources Association, v. 41, 

no. 3, p. 621–633.

Maechler, Martin; Rousseeuw, Peter; Struyl, Anja; Hubert, 

Mia; and Hornik, Kurt, 2015, Cluster—Cluster anaysis 

basics and extensions: R Foundation for Statistical Computing, R package version 1.15.2

McLeod, A.I., 1994, Diagnostic checking of periodic 

autoregression models with application: Journal of 

Time Series Analysis, v. 15, p. 221–233, accessed 

August 20, 2015, at http://onlinelibrary.wiley.com/

doi/10.1111/j.1467-9892.1994.tb00186.x/pdf.

Menne, M.J., Williams, C.N., Jr., and Vose, R.S., 2014, United 

States Historical Climatology Network (USHCN) Version 

2.5 serial monthly dataset: National Climatic Data Center, National Oceanic and Atmospheric Administration, 

accessed May 15, 2014, at http://cdiac.ornl.gov/ftp/ushcn_

v2.5_monthly/.

Milly, P.C.D., Dunne, K.A., and Vecchia, A.V., 2005, Global 

pattern of trends in streamflow and water availability in a 

changing climate: Nature, v. 438, no. 7066, p. 347–50.

Murphy, E.C., Fritz, A.K., and Fleming, R.F., 1997, The 

Jerusalem and Tolna outlets in the Devils Lake Basin, North 

Dakota: North Dakota Geological Survey Report of Investigation No. 100, 36 p.
48 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

National Oceanic and Atmospheric Administration, National 

Climatic Data Center, 2014, Tree ring: National Oceanic 

and Atmospheric Administration, accessed July 15, 2014, 

at http://www.ncdc.noaa.gov/data-access/paleoclimatologydata/datasets/tree-ring.

Pettyjohn, W.A., 1967, Geohydrology of the Souris River Valley in the vicinity of Minot, North Dakota: U.S. Geological 

Survey Water-Supply Paper 1844, 53 p.

Pusc, S.W., 1994, Hydrogeology of the Minot Aquifer, Ward 

County, North Dakota: North Dakota State Water Commission Ground-water Studies Number 102, Part 2, 92 p.

R Core Team, 2015, R—A language and environment for 

statistical computing: Vienna, Austria, R Foundation for 

Statistical Computing, September 30, 2014, at http://

www.R-project.org.

Rannie, W.F., 1998, A survey of hydroclimate, flooding, and 

runoff in the Red River Basin prior to 1870: Geological 

Survey of Canada Open-File Report 3705, 189 p.

Rao, L.Y., Sun, Ge, Ford, C.R., and Vose, J.M., 2011, Modeling potential evapotranspiration of two forested watersheds 

in the southern Appalachians: Transactions of the ASABE, 

v. 54, no. 6, p. 2067–2078.

Rosenberry, D.O., Stannard, D.I., Winter, T.C., and Martinez, 

M.L., 2004, Comparison of 13 equations for determining evapotranspiration from prairie wetland, Cottonwood 

Lake area, North Dakota, USA: Wetlands, v. 24, no. 3, p. 

483–497.

Ryberg, K.R., 2015, The impact of climate variability on 

streamflow and water quality in the North Central United 

States: Fargo, North Dakota State University, Ph.D. dissertation, 277 p.

Ryberg, K.R., Lin, Wei, and Vecchia, A.V., 2014, Impact of 

climate variability on runoff in the north-central United 

States: Journal of Hydrologic Engineering, v. 19, no. 1, p. 

148–158.

Shapley, M.D., Johnson, W.C., Engstrom, D.R., and Osterkamp, W.R., 2005, Late-Holocene flooding and drought in 

the Northern Great Plains, USA, reconstructed from tree 

rings, lake sediments, and ancient shorelines: The Holocene, 

v. 15, p. 29–41.

Small, David, and Islam, Shafiqul, 2008, Low frequency variability in fall precipitation across the United States: Water 

Resources Research, v. 44, no. 4, 13 p.

Small, David, and Islam, Shafiqul, 2009, A synoptic view of 

trends and decadal variations in autumn precipitation across 

the United States from 1948 to 2004: Journal of Geophysical Research, v. 114, no. D10, 14 p.

St. George, S., and Nielsen, E., 2002, Hydroclimatic change in 

southern Manitoba since A.D. 1409 inferred from tree rings: 

Quaternary Research, v. 58, no. 4, p. 103–111.

St. George, S., and Nielsen, E., 2003, Palaeoflood records for 

the Red River, Manitoba, Canada, derived from anatomical 

tree-ring signatures: The Holocene, v. 13, no. 4, p. 547–555.

U.S. Army Corps of Engineers, 2013, Regional and reconstructed hydrology of the Souris River: St. Paul District, 

158 p.

U.S. Department of Agriculture, 2014, Gridded Soil Survey Geographic Database (gSSURGO) Version 1.1: 

accessed October 21, 2014, at https://gdg.sc.egov.usda.gov/

GDGOrder.aspx.

U.S. Department of Commerce, 2012, The Missouri/Souris 

floods of May–August 2011: Kansas City, National Weather 

Service, 68 p.

U.S. Geological Survey, 2014a, The National Map Viewer: 

accessed September 16, 2014, at http://viewer.nationalmap.

gov/viewer/.

U.S. Geological Survey, 2014b, USGS water data for North 

Dakota: U.S. Geological Survey, National Water Informatoin System Web page, accessed November 18, 2014, at 

http://waterdata.usgs.gov/nd/nwis/. 

U.S. Geological Survey, 2014c, StreamStats Version 3: 

accessed October 8, 2014, at http://ssdev.cr.usgs.gov/

ss_dev/.

U.S. Geological Survey, 2015, National Hydro Dataset 

(NHD0901): U.S. Geological Survey, The National Map, 

accessed September 8, 2014, at http://nhd.usgs.gov/wbd.

html.

Vecchia, A.V., 2002, Simulation of a proposed emergency 

outlet from Devils Lake, North Dakota: U.S. Geological 

Survey Water-Resources Investigations Report 02–4042, 

129 p.

Vecchia, A.V., 2008, Climate simulation and flood risk 

analysis for 2008–40 for Devils Lake, North Dakota: U.S. 

Geological Survey Scientific Investigations Report 2008–

5011, 28 p., accessed June 23, 2015, at http://pubs.usgs.gov/

sir/2008/5011/.

Wiche, G.J., Vecchia, A.V., Osborne, L., Wood, C.M., and Fey, 

J.T., 2000, Climatology, hydrology, and simulation of an 

emergency outlet, Devils Lake Basin, North Dakota: U.S. 

Geological Survey Water-Resources Investigations Report 

00–4174, 16 p.
References Cited 49

Whitcher, Brandon, 2013, Waveslim—Basic wavelet routines 

for one-, two- and three-dimensional signal processing: R 

Foundation for Statistical Computing, R package version 

1.7.3, accessed October 6, 2014, at http://CRAN.R-project.

org/package=waveslim.

Yang, Song, Ding, X., Zheng, D., and Li, Q., 2007, Depiction 

of the variations of Great Plains precipitation and its rela-

tionship with Tropical Central-Eastern Pacific SST: Journal 

of Applied Meteorology and Climatology, v. 46, no. 2, p. 

136–153.

Yonetani, T., and Gordon, H.B., 2001, Abrupt changes as 

indicators of decadal climate variability: Climate Dynamics, 

v. 17, no. 4, p. 249–58.

Appendix 51

Appendix. Water-Balance Model Equations

The water-balance model (WBM) consists of a series of equations, which, when operated in sequential order, estimate the 

amount of available water, either stored in the soil profile or contributing to runoff, at a monthly time step. Model inputs are indepth equivalents and consist of monthly precipitation, average temperature, and potential evapotranspiration (PET) along with 

static variables such as available water storage capacity (AWSC) and soil permeability (Ks). Model outputs estimate the amount 

of available water contributing to runoff and are converted from depth equivalent values to streamflow estimates in cubic meters 

per second. 

The first parameters in the model are two coefficients for adjusting PET computed using the Hamon method, as described 

in the section “Water-Balance Model Description and Calibration”: 

PET

C PET

C PET

PET

M H

J H

H

=





 1 1

1 1

1 1

. ,

. ,

. ,

May

June

remaining months 













where 

PET is adjusted PET (millimeters), 

PETH is PET computed using the Hamon method (millimeters), 

CM is the adjustment factor PET (May) found in table 8 (dimensionless), and 

CJ is the adjustment factor PET(June) found in table 8 (dimensionless). 

The next model parameter is a coefficient for adjusting the available soil-moisture storage capacity computed using the 

geographic information system (GIS) data described in the “Water-Balance Model Description and Calibration” section of this 

report:

AWSC = C A AWS G WSC IS

where

AWSC is the adjusted available water storage capacity (millimeters), 

AWSCGIS is the storage capacity computed from the GIS data (millimeters), and 

CAWS is the available water storage adjustment coefficient (dimensionless). 

As indicated in table 8, for most of the subbasins this coefficient (CAWS ) was equal to one; however, some of the coefficients 

were larger than one, signifying depression storage may play a larger role in the amount of water retained within the subbasin 

as compared to running off, and some of the coefficients were less than one, signifying a subbasin with less than normal depression storage. One subbasin, Souris River near Verendrye, North Dakota (05120000), had a coefficient of 0.9, which may be the 

consequence of more runoff (less infiltration) being delivered from urban/suburban development in this particular subbasin that 

contains the city of Minot.

The other WBM parameters and equations depended on how precipitation entered the system (rain or snow). The classification of precipitation as rain or snow is dependent on monthly average temperature (Tavg) such that precipitation is considered rain 

if Tavg exceeds a specified threshold (Tr), snow if Tavg is less than a specified threshold (Ts), and a mixture of rain and snow if Tavg

is between the two thresholds: 

P

P T T

P

T T

T T

T T T

T T

snow

avg s

r avg

r s

s avg r

avg r

=

≤ = −

−

− < <

≥ =

,

,

,

10

0 2

 Co

o

C




























52 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

P P rain = − Psnow

where

Psnow is snow (millimeters), 

P is monthly precipitation (millimeters), and 

Prain is rain (millimeters). 

The values selected for Tr (2 degrees Celsius [oC]) and Ts (-10 oC) were similar to those selected by Gray and McCabe (2010), 

who used 5 oC for Tr and minus 5 oC for Ts. Differences in the thresholds for different climatic regions, topography, vegetation, 

and so forth, would be expected for a monthly time step. Gray and McCabe (2010) studied the Upper Yellowstone River in the 

Rocky Mountains. 

Precipitation classified as snow is accumulated in snowpack (SNp). Subtractions from SNp only take place when Tavg is 

greater than (>) Ts and are lost in the form of snowmelt (SM). The fraction of SNp lost to SM is directly related to the difference 

between Tavg and Ts as well as the time of year:

SM

T T

T T SN T T

avg s

= avg s P avg s

≤

− > −

0

10

,

min[ ( ), ], (for January March)

min[20( ) T T , ] SN ,T T (for April December) avg s

− > P avg s −





















where

SM is the snowmelt (millimeters) and

SNp is the accumulated snowpack (millimeters). 

The smaller value of the multiplier (10) for January–March was selected through trial and error to match the timing of the 

snowmelt season and is reasonable considering more energy (higher temperature) would be required to condense and ripen the 

snowpack early in the snowmelt season than during later months. A small portion of SM is assumed to become direct runoff during the late winter and early spring months when the ground is still frozen/thawing and not all of the SM infiltrates the soil:

SM

C SM

SM C SM

RO

SM

SM =

−

+

min( , ),

min( . , ),

1

0 01 5

 January February

 March

 April

 May December

min( , ),

,

C S SM M 15

0 −





























where

SMRO is the direct snowmelt runoff (millimeters) and 

CSM is the snowmelt adjustment coefficient (dimensionless; table 8). 

A SMRO cap of 1 millimeter (mm) was used between January and February, 5 mm during March, 15 mm during April, and no 

SMRO was considered between May and December (SMRO = 0). A cap was placed on SMRO such that monthly values would not 

exceed reasonable values for the region.

In addition to SMRO, groundwater runoff also contributes directly to streamflow, but unlike SMRO, groundwater runoff is 

subtracted from the initial soil moisture remaining from the end of previous time step:

GW F AWS RO = GW i

AWS A u i = − WS GWRO

where

GWRO is groundwater runoff (millimeters), 

FGW is the fraction of soil moisture that becomes groundwater runoff (dimensionless), 

AWSi is initial available soil moisture from the previous time step (millimeters), and 

AWSu is updated (for the current time step) available soil moisture (millimeters). 
Appendix 53

FGW depends on temperature, permeability, and the degree of soil moisture deficit:

F

T T T

T T

AWS

AWSC

e GW

avg r s

r s

i = −

−





 









 



 0 02 0

2 1 4

. max

min( , ) , . min(K , ) s 20

20 −1 



 





where

Ks is the soil permeability (centimeters per hour). 

Ks was computed using the GIS data described in the “Water-Balance Model Description and Calibration” section. The constant 

multiplier (0.02) was determined through the calibration process. The second term on the right-hand side equals one if Tavg is 

greater than or equal to (≥) Tr, zero if Tavg is less than or equal to (≤) Tr , and lies in between 0 and 1 otherwise, and is present to 

reduce the groundwater runoff during winter months when soils are frozen. The third term equals one if AWSi equals AWSC and 

decreases as AWSi approaches zero. The last term on the right-hand side is one if Ks≥20 centimeters per hour and decreases to a 

small value as Ks tends towards zero. 

For the months of April and May, when there is both ground thaw and snowmelt, GWRO determined from the previous equations is increased up to a maximum potential value and the extra groundwater runoff is taken from snowmelt: 

GW

GW SM SM GW GW

GW SM SM SM SM GW G RO

RO RO

RO RO RO

* max max

max

,

, = − ≥ −+ − − < − WRO













 for April − May

SM SM SMRO GWRO GWRO

* * = − − − ( ) for April − May

GWRO GWRO

* = − , for January March and June − December

SM SM SMRO

* = − , for January − − March and June December

where

GWRO

* is the adjusted groundwater runoff (millimeters),

GWmax is determined by replacing AWSi with AWSC in the previous equation for FGW and multiplying the resulting 

fraction times AWSC, and 

SM *

 is the current snowmelt remaining after subtracting snowmelt runoff and (or) extra groundwater runoff 

(millimeters). 

After snowmelt contributions to groundwater have been determined, the remaining SM * is combined with the remaining 

soil-moisture storage and rainfall input to obtain the total available water in soil-moisture storage:

AWS A u u WS SM Prain

* * = + +

where

AWSu

* is the current total available water in soil-moisture storage (millimeters). 

If AWSu

* exceeds AWSC, more direct runoff is taken from the surplus before considering water losses to actual evapotranspiration (AET): 

Sur1 = − max(AWS A u WSC, ) 0 *
54 Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow

DS

C Sur

C Sur

C Sur RO

DRO

DRO

DRO

=

0 65

0 05

0 01

1

1

1

.

.

.

 , April

 , May

 , October November

 otherwise

−



























0, 

Sur Sur DS 2 1 = − RO

where

Sur1 is the initial surplus (millimeters), 

DSRO is the direct runoff from the initial surplus (millimeters), 

CDRO is the direct runoff adjustment coefficient (dimensionless; table 8), and 

Sur2 is the surplus remaining after direct runoff (millimeters). 

The direct runoff coefficient varies by month such that April has the largest coefficient, whereas May, October, and November 

have much smaller coefficients. In April most of the ET comes later in the month, when surplus would be expected to runoff 

before ET has a chance to begin. Similarly, in May, there tends to be somewhat higher ET late in the month compared to early 

in the month. In October and November the ground is beginning to freeze and the crops are becoming dormant, which limits the 

amount of water retained in the soil profile and results in slightly more runoff. 

The remaining surplus (Sur2), if any, is used to satisfy primary ET (the amount of ET that does not come from the soilmoisture storage) and if primary ET does not satisfy PET, additional ET is removed from the soil-moisture storage: 

AET

PET Sur PET

PET Sur F AWSC Sur PET

PET F

ST

S

=

≥

+ < <

,

min( , ),

min( ,

2

2 2 0

T u AWS Sur *

), 2 ≤ 0





















F T ST = + avg

1

25

max(min( , ) 24 1 0, )

where 

AET is actual ET (millimeters) and 

FST is the fraction of soil-moisture storage that is exposed to secondary ET (dimensionless). 

This fraction is dependent on temperature, such that temperatures below -1oC result in none of the soil moisture being available 

for secondary ET and as temperatures increase so does the amount of soil-moisture storage available for secondary ET. If Tavg

exceeds 24 oC, all of the soil moisture can potentially be lost to secondary ET. The reason for the temperature dependence is 

that most of the secondary ET is considered to result from evapotranspiration from natural vegetation and crops, which would 

become dormant as temperatures near freezing and grow rigorously during mid- to late summer. 

After water losses to AET are considered, remaining water surplus (if any) goes to excess overland runoff:

EORO = − max(Sur AET, ) 2 0

where

EORO is excess overland runoff for the current time step (millimeters). 

To account for time to travel across the watershed, a lag was introduced such that a portion of the excess overland runoff was 

assumed to enter the stream in the current month and a portion in the subsequent month (table 8):

EORO F E EO i R O F O i EO i R EO O i

* = + , , − − 1 1 , ,
Appendix 55

where

EORO

* is the current excess overland runoff (millimeters), 

FEO,i−1 is the fraction of excess overland runoff from the previous month entering the stream for the current month 

(dimensionless), 

EORO,i−1 is excess overland runoff (millimeters) from the previous month,

FEO,i is the fraction of excess overland runoff for the current month entering the stream (dimensionless), and 

EORO,i is the excess overland runoff (millimeters) for the current month. 

Considering all of the previous water gains and losses in aggregate, the total runoff for the current month is

TOT S RO = + M G RO W D RO + + S E RO ORO

* *

where 

TOTRO is the total runoff (millimeters).

The remaining storage at the end of the current time step (and the beginning storage at the start of the next time step) is

AWS A i u +1 = − WS AET − − DSRO EORO i

*

,

 where

AWSi+1 is the remaining storage at the end of the current time step (millimeters).

Once a depth equivalent of TOTRO is determined for each grid cell, an average is taken over the entire subbasin and a streamflow equivalent calculated:

Q Avg TOT A d = WS RO WS 1, { 000 } / ( , 86 400 )

where

Q is monthly mean streamflow (cubic meters per second); 

AvgWS{.} is the average of the total runoff for the grid cells in the watershed (millimeters); 

AWS is the watershed area (square kilometers); and 

d is the number of days in the month.

The equation has a unit correction factor of 1,000 to convert runoff from a depth equivalent to volume in cubic meters and a unit 

correction factor of 86,400 to convert from days to seconds.
Publishing support provided by: 

Rolla Publishing Service Center 

For more information concerning this publication, contact:

Director, USGS North Dakota Water Science Center

821 East Interstate Avenue

Bismarck, North Dakota 58503

(701) 250–7400

Or visit the North Dakota Water Science Center Web site at: 

http://nd.water.usgs.gov/

Kolars and others—Stochastic Model for Simulating Souris River Basin Precipitation, Evapotranspiration, and Natural Streamflow—SIR 2015–5185

ISSN 2328-0328 (online)

http://dx.doi.org/10.3133/sir20155185