do not necessarily reflect the views of UKDiss.com.

**Evaluation and Mitigation of the Temporal Evolution of Microbial Contamination Risk in Surface Water Systems**

**Abstract**

This thesis presents a methodology for evaluating the temporal evolution of microbial contamination risk in surface water systems. The objective of this research is to advance the field of risk assessment across the realm of environmental engineering by creating a methodology that is applicable to numerous complex environmental phenomena, as well as developing the tools and theories necessary to make this methodology pragmatically viable. This is accomplished in part by developing a probabilistic mass conservation model to describe large scale dynamics of microbial contaminants from nonpoint sources throughout a watershed, in a manner that embraces the uncertainty inherent to these phenomena. With a developed microbial fate and transport model that attributes contamination to discrete spatial sources and pathways, the action that must be taken to minimize this contamination risk can be determined. Any method utilized to minimize nonpoint source microbial contamination risk will have an associated cost of implementation, so an optimization algorithm is applied to the results of the probabilistic mass conservation model to determine the optimal placement of control measures throughout the watershed necessary to achieve an objective reduction in contamination risk. To aid in the practicality of risk assessment within complex environmental systems, a rapidly converging algorithm is developed to efficiently propagate uncertainty through complex systems, such as that which is described by the probabilistic mass conservation model. This describes the time and space evolving probability density function of microbial contamination, which in turn is used to determine how risk changes throughout time. Finally, a standardized methodology for parameterizing the probabilistic mass conservation model is presented, which in turn allows for the derivation of certain transition rates using statistical and physical first principles. The ability of environmental engineers and watershed managers to perform a large-scale, systems-driven assessment of risk in environmental systems has not been previously viable, due to the irreducible complexity of environmental processes and the significant computational cost demanded by conventional means of uncertainty propagation for risk assessment. This project provides a versatile, scalable platform for the assessment of the temporal evolution of risk in environmental systems, applied specifically to microbial contamination in surface water systems.

**Preface**

In Chapter 1, a methodology is presented to determine the temporal evolution of microbial contamination risk in surface water systems. To describe the temporal evolution of risk of microbial contamination, system components must first be identified. System components are the potential sources and pathways of microbial contamination that contribute to water quality failure. These components are identified by first relating each spatial segment of the watershed to its contribution of microbial contaminants to the surface water system using sets of directed paths. This captures the influence that spatially dependent characteristics of the watershed have on nonpoint source microbial contamination. Next, the stochastic spatial-temporal behavior of microbial contamination is described within each component using a probabilistic mass conservation model. Finally, the system components must be recombined within the system reliability framework to create the metric of temporal evolution of risk. Using the temporal evolution of risk, sustainability that aid in the interpretation of the temporal evolution of risk metrics will be created.

In Chapter 2, a linear optimization algorithm is applied to results of probabilistic mass conservation model to determine the optimal implementation of control measures and the reduction of contamination risk. This can be summarized in three different schemas: single point control, multiple point control, and multiple point control through time.

In Chapter 3, an algorithm is developed to efficiently propagate uncertainty through complex systems, such as that which is described by the probabilistic mass conservation model. The algorithm allows the description of time and space evolving probability density function of microbial contamination, which in turn is used to determine how risk changes throughout time. The general function of this algorithm is to determine an output uncertainty by identifying the combination of input values that result in an output that is one standard deviation above the mean output, where the output standard deviation is unknown. This is accomplished by numerically approximating the influence that each parameter has on the output, then using this influence in tandem with a converging algorithm to determine the standard deviation of the output.

In Chapter 4, a methodology is presented for parameterizing the probabilistic mass conservation model, specifically through deriving transition rates using statistical and physical first principles. The importance of treating transition rates in this manner is demonstrated by showing that a number of transition rates can be derived using statistical first principles to determine the conditional transition probability. It is shown that while such derivations provide elegant formulas for certain transition rates, they can become extremely complex for other parameters. However, this increase in complexity does not preclude the importance of treating transition rates as conditional probabilities divided by a time increment. The Bayesian nature of these transition rates has a bearing on the empirical approximation of transition rates that have unwieldly purely theoretical derivations, and provides insight into the probabilistic distribution of these transition rates.

** **

This page Left Intentionally Blank

Chapter 1: Temporal Evolution of Risk and Sustainability Metrics

1.2.1 Reachability and Directed Paths

1.2.2 Microbial Transport along Directed Paths

1.2.4 Temporal Evolution of Risk

1.2.5 Sustainability as a Function of Time

1.3.1 Temporal Evolution of Risk Results

1.3.4 Sustainability Index Results

1.4.1 Combination of Metrics for Engineering Judgement

1.4.2 Use of Metrics in Decision Making

1.4.3 Potential Applications in System Controllability

Chapter 2: Contributor Based Control of Microbial Contamination

2.2 Proportional Reduction in Contributor Based Control

2.2.1 Proportional Reduction in Single Point Control Schema

2.2.2 Single Point Control Application

2.2.3 Multiple Point Control Schema

2.2.4 Multiple Point Control Application

2.2.5 Multiple Point Control through Time Schema

2.3 Control Application Summary

2.3.2 Possible Means of Control through Application of Management Practices

2.3.3 Grass Buffer Strips Effectiveness

2.3.4 Grass Buffer Strips Feasibility

2.3.5 Source Control Effectiveness

2.4.1 Single Point Control Results

2.4.2 Single Point Control Discussion

2.4.3 Multiple Point Control Results

2.4.4 Multiple Point Control Discussion

2.4.5 Multiple Point Control through Time Results

2.4.6 Multiple Point Control through Time Discussion

Chapter 3: Rapidly Converging Uncertainty Propagation Algorithm

3.2 Uncertainty Propagation Theory

3.4 Uncertainty Propagation with Correlated Random Variable Inputs

Chapter 4: Parameterization of Probabilistic Mass Conservation Model

4.2 Probabilistic Nature of Transition Rates

4.3 Derivation of Transition Rates

4.3.1 Lateral Transport Transition Rate Substantiation

4.3.2 Settling Transition Rate Substantiation

4.3.3 Uptake Transition Rate Substantiation

4.3.4 Empirical Derivation of Parameters Methodology

4.4 Probabilistic Distribution of Transition Rates

Appendix A: Probabilistic mass conservation model MATLAB Code

Appendix B: Hill Slope and Stream Tables

Appendix C: Temporal Evolution of Risk MATLAB Code

Appendix D: Control Optimization MATLAB Code

Appendix E: Rapid Uncertainty Propagation Algorithm MATLAB Code

Appendix F: Temporal Sustainability Index MATLAB Code

**Index of Figures**

** **

Figure 1.1.A: Example of Simple Watershed

Figure 1.1.B: Extracted Dendritic Network

Figure 1.2: Schematic of Resulting System Components

Figure 1.3: Microbial Transport Model

Figure 1.4: Microbial Transport Finite Volume Method

Figure 1.5: Fault Tree Diagram

Figure 1.6: Collapsed Fault Tree Diagram

Figure 1.7.A: Shepherd creek watershed

Figure 1.7.B: Sub-Catchment of Shepherd Creek Watershed

Figure 1.8.A: Concentration Mean Curve of System Components

Figure 1.8.B: Concentration Standard Deviation Curve of System Components

Figure 1.9: Concentration Mean Curve of System with Measured Data

Figure 1.11: Conventional Resilience

Figure 1.13: Vulnerability Results

Figure 1.14: Sustainability Index Results

Figure 2.1: Sample Watershed with Objective Points

Figure 2.2: Concentration Mean Curve at Point B

Figure 2.3: Concentration Mean Curve at Point C

Figure 3.1: Convergence of Guessed and Actual Standard Deviation

Figure 3.2: Convergence of Uncertainty Propagation Algorithm

Figure 3.3: Comparison of Uncertainty Propagation Algorithm and MCS for Time Dependent System

Figure 3.4: Mean Solution of Lorenz Equation

Figure 3.5: Comparison of Uncertainty Propagation Algorithm and MCS for Time Dependent System

Figure 4.1: Free Body Diagram of Soil Particle on Hill Slope

Figure 4.2.A: Sampling Results from Beta Distribution with

clear x

for it=1:length(std_array_new)

%Evaluate the Mean and standard deviation from each component at each point in time

mean_at_time=mean_array_new(it,:);

std_at_time=std_array_new(it,:);

%Evaluate the combined mean and standard deviation from all components at

%each point in time

m=sum(mean_at_time);

s=sqrt(sum(std_at_time.^2));

%This improves stability for very low probability of failures

if s<1e-11

Prob(it)=0;

else

%Determine parameters of lognormal distribution

mu = log((m^2)/sqrt(s^2+m^2));

sigma = sqrt(log(s^2/(m^2)+1));

%Determine probability of failure through cumulative lognormal

%distribution.

Prob(it)=1-logncdf(Limit,mu,sigma);

end

end

**Appendix D: Control Optimization MATLAB Code**

%Controlability

clear all

clc

N=11; %Number of Nodes

C=9; %Number of Control Points

%Import Contamination Data for Control Points based on Probabilistic Mass

%Conservation Model

L=transpose(xlsread(‘Shep Control.xlsx’,’Sheet3′,’B8:J18′));

%This is our goal load at the point of interest, so the accumulation of all

%the individual nodes should be below this.

Lg=5e4;

%Proportional reduction

%So this portion creates the matrix we are going to use. The decision

%variable X is the proportion that we reduce each node load by.

for ic=1:C

for i=1:N

A(ic,i)=-L(ic,i);

end

end

%Establish Initial Constraint Matricies

AB=zeros([2*N+C N]);

AB(C+1:N+C,:)=eye(N,N);

AB(N+C+1:2*N+C,:)=-eye(N,N);

%Establish constraint that sum of the loads minus their reduction to be less than the

%goal.

for ic=1:C

AB(ic,1:N)=A(ic,1:N);

B(ic)=Lg-sum(L(ic,:));

end

%Determine Maximum proportional reduction

B(C+1:N+C)=0.65;

B(N+C+1:2*N+C)=0;

B=transpose(B);

z(1:N)=1;

%Perfom Linear Optimization to determine optimal proportional reduction

[X Z]=linprog(z,AB,B)

**Appendix E: Rapid Uncertainty Propagation Algorithm MATLAB Code**

%Rapid Uncertainty Propagation Algorithm

clear all

%Establish Inputs

s=rng(3);

Vars=4;

InputMean=rand(Vars,1)*1;

Inputstd=rand(Vars,1)*0.1;

%Initialize count

Count=1;

%Specify Time of Evalutation (if relevant)

eval_time=500;

%Create Vector of Input Means

In=InputMean;

%Run function with output gC

GcFun

Count=Count+1;

%Because the mean inputs were entered into the model, the output is the

%mean output

Mean=gC;

%Specify Variable Number

var_num=length(InputMean);

%Create range of input weights from -3 to 3, create matrix were each weight

%is varied across this range while all others are 0.

gindex=1;

range=3;

clear w

for q=1:var_num

for iw=-range:range

w(gindex,1:2)=zeros(1,2);

w(gindex,q)=iw;

gindex=gindex+1;

end

end

%Sweep across the input weight matrix

for iW=1:length(w)

%Apply input weights to inputs

In=InputMean+Inputstd.*transpose(w(iW,:));

%Run Model

GcFun

Count=Count+1;

%Save Results

DiffOutputs(iW)=gC;

InSave1(:,iW)=In;

end

y=DiffOutputs;

%Find vector of differences between output and mean

D=(DiffOutputs-Mean);

%For stability

for itna=1:length(D)

if isnan(D(itna))==1

D(itna)=0;

end

end

%Determine Response Space Characteristics

syms x

iRc=1;

%Determine Dimension Size

Dimsize=length(w)-6;

for iw=1:7:Dimsize

%Perform linear regression on response space

[polyvec]=polyfit(InSave1(iRc,iw:iw+6)-InputMean(iRc),D(iw:iw+6),1);

%Save first parameter of linear regression

Mvec(iRc)=polyvec(1);

%Determine correlation coefficient between realizations across

%characteristic of response space

M=corrcoef(-3:1:3,D(iw:iw+6));

%Compare to correlation coefficient between differences across

%characteristic of response space

MP=corrcoef(InSave1(iRc,iw:iw+6)-InputMean(iRc),D(iw:iw+6));

%Extract Correclation Coefficient

Rc(iRc)=M(2,1);

RcP(iRc)=MP(2,1);

%Determine the influence vector

Infl(iRc)=mean(diff(D(iw:iw+6)));

%Indexing

iRc=iRc+1;

end

%Allow for iteration of guesses

for iCC=1:1

%Save initial Influence guess

Infl_f=Infl;

%Solve for adjusted influence

iW=1;

%Clear as necessary

clear w

clear A DiffOutputs StdGuess

DiffOutputs(1)=gC;

%Create initial Guess

StdGuess(1)=Mean;

%Save Guess

Std_Record(1)=StdGuess(1);

Esave=1;

%Determine convergence threshold

Thresh=Mean*1e-04;

%While system is above this threshold…

while abs(DiffOutputs(iW)-Mean-Std_Record(iW))>Thresh

%Based on guessed influence and guessed STD, generate vector of guessed

%input weights

w=Infl./StdGuess(iW);

%Apply guessed input weight vector to inputs

In=InputMean+Inputstd.*(w(:));

%Run Model

GcFun

%Save Output

Count=Count+1;

DiffOutputs(iW+1)=gC;

%Evaluate Local response space gradients about the guessed input weights

iRc=1;

for iw=1:7:Dimsize

%Spline interpolate locally

NewD(iw:iw+6)=interp1(-3:3,D(iw:iw+6),[w(iRc)+3*w(iRc)/-10]:w(iRc)/10:w(iRc)+3*w(iRc)/10,’spline’);

Infls(iRc)=mean(diff(NewD(iw:iw+6)))/w(iRc)*10; %Local RS gradient

iRc=iRc+1;

end

%Record Outputs

Std_Record(iW+1)=StdGuess(iW);

w_record(:,iW+1)=w;

iW=iW+1;

%See difference between initial guessed influence and local response space

%gradient

E1=Infls-Infl_f;

%Adjust Influence guess by error

Infl=Infl-1/20*E1;

%Save cumulative error

Esave(iW)=sum(E1.^2);

%Backpropagate error to determine new standard deviation guess.

StdGuess(iW)=abs(StdGuess(iW-1)-.1*(StdGuess(iW-1)-(DiffOutputs(iW)-Mean)));

end

%Save Results

True_STD1=Std_Record(iW);

w_true_z=Infl./StdGuess(iW);

Infl_z=Infl;

Efinal(iCC)=Esave(iW);

%Adjust new guessed influence by stochastic gradient descent of error

%between the influence that produces weights with the correct local gradient

%and the guessed influence

Infl=Infl_f-(Infl_f-Infl_z).*(1-Rc.^2)

%Determine Standard Deviation For fixed, new influence guess

iW=1;

%Clear as necessary

clear w

clear A DiffOutputs StdGuess

%Initialize with new guesses.

DiffOutputs(1)=gC;

StdGuess(1)=Mean;

Std_Record(1)=StdGuess(1);

%While system is above this threshold…

while abs(DiffOutputs(iW)-Mean-Std_Record(iW))>Thresh

%Based on guessed influence and guessed STD, generate vector of guessed

%input weights

w=Infl./StdGuess(iW);

%Apply guessed input weight vector to inputs

In=InputMean+Inputstd.*(w(:));

%Run Model

GcFun

%Save Outputs

Count=Count+1;

DiffOutputs(iW+1)=gC;

%Save Outputs

Std_Record(iW+1)=StdGuess(iW);

w_record(:,iW+1)=w;

iW=iW+1;

%Save Error of this guess

Esave(iW)=sum(E1.^2);

%Backpropagate Error to determine new standard deviation guess

StdGuess(iW)=abs(StdGuess(iW-1)-.1*(StdGuess(iW-1)-(DiffOutputs(iW)-Mean)));

end

%We’ve converged on the correct standard deviation

True_STD(iCC)=StdGuess(iW);

end

**Appendix F: Temporal Sustainability Index MATLAB Code**

%Temporal Sustainability Metrics

clc

clear all

%Import Results of Uncertainty Propagation Algorithm for Microbial

%Contamination

stdarray=xlsread(‘Array.xlsx’,’Sheet1′);

meanarray=xlsread(‘Array.xlsx’,’Sheet2′);

%Assign Safety Limit

LimitB=6.2e4;

iL=1

%Create Vector of various different Limits

LimitVector=(LimitB-LimitB/2:LimitB/15:LimitB+LimitB/2)/283.168;

%Evaluate metrics for different Limits

for Limit=LimitB-LimitB/2:LimitB/15:LimitB+LimitB/2;

%Evaluate Temporal Evolution of Risk

for it=1:length(stdarray)

%Extract the mean and standard deviation of each component at this

%point in time

mean_at_time=meanarray(it,:);

std_at_time=stdarray(it,:);

%Combine mean and standard deviation for all system components

m=sum(mean_at_time);

s=sqrt(sum(std_at_time.^2));

%For stabilitiy.

if s<0

Prob(it)=0;

else

%Determine parameters of lognormal distribution

mu = log((m^2)/sqrt(s^2+m^2));

sigma = sqrt(log(s^2/(m^2)+1));

%Determine probability of failure through cumulative lognormal

%distribution.

Prob(it)=1-logncdf(Limit,mu,sigma);

end

%Save Results

mean_sum(it)=m;

std_sum(it)=s;

end

iR=1;

%Determine for different Resilience Limits of acceptable risk

for Res_Limit=.01:.05:1;

%Vectorize Resilience Limit

ResLine(1:length(Prob))=Res_Limit;

%Find Intersect of temporal evolution of risk curve with acceptable

%risk line

IntersectX=InterX([1:length(Prob);ResLine],[1:length(Prob);Prob]);

%For Stability, make zero if empty

if isempty(IntersectX)==1

Resilience(iL,iR)=0;

Integral(iL,iR)=0;

else

%Clear as necessary

clear durations

inX=1;

%Determine duration of unacceptable probability of failure

for interii=2:2:length(IntersectX)

%An exceedance will need to have two intersects

durations(inX)=IntersectX(1,interii)-IntersectX(1,interii-1);

inX=inX+1;

end

%Determine Intersection area above curve

IntersectionsRes(iL,iR,:)=[int32(IntersectX(1,1)),int32(IntersectX(1,length(IntersectX)))];

%Determine resilience for this acceptable risk and this contamination limit

Resilience(iL,iR)=sum(durations);

%Check the Integral

Integral(iL,iR)=trapz(Prob(int32(IntersectX(1,1)):int32(IntersectX(1,length(IntersectX))))./100);

end

iR=iR+1;

end

%Vulnerability is the maximum probability of failure in respect to a

%certain threshold

Vulnerability(iL)=max(Prob);

%Save Probability

ProbMatrix(:,iL)=Prob;

iL=iL+1

end

iL=2;

iR=2;

%Convert Resilience and Vulnerability into format that can be used to

%determine sustainability index

for it=1:length(ProbMatrix(:,iL))

if it<=IntersectionsRes(iL,iR,1)

Res_time(it)=Resilience(iL,iR);

elseif it>IntersectionsRes(iL,iR,1) && it<=IntersectionsRes(iL,iR,2)

Res_time(it)=Resilience(iL,iR)-(it-IntersectionsRes(iL,iR,1));

else

Res_time(it)=0;

end

Vul_time(it)=max(ProbMatrix(1:it,iL));

end

%Determine sustainability index

SI=(1-ProbMatrix(:,iL)).*(max(Res_time)-Res_time(:))/100.*(1-Vul_time(:));

%Adjust sustainability index

SI_true=transpose(1./(.01.*Resilience(:,3))).*(1-Vulnerability)*(1-Prob(10));

hold on

%Plot results in CFU/100mL

plot((1:length(stdarray))/100,stdarray/283.168)

plot((1:length(std_sum))/100,std_sum/283.168,’k’)

xlabel(‘Time (Hour)’)

ylabel(‘Concentration (CFU/100mL)’)

title(‘Concentration Standard Deviation Curve (7/11/06)’)

figure

hold on

plot((1:length(meanarray))/100,meanarray/283.168)

plot((1:length(mean_sum))/100,mean_sum/283.168,’k’)

xlabel(‘Time (Hour)’)

ylabel(‘Concentration (CFU/100mL)’)

title(‘Concentration Mean Curve (7/11/06)’)

data=[0 0 4500 13500 60000 50000 20000 0];

time=[’12:00′ ’15:00′ ’16:00′ ’19:00′ ‘1:00’ ‘4:00’ ‘7:00′ ’10:00’];

hour=[1 4 5 8 14 17 20 23];

data_time=hour/.01;

figure

plot(.01:.05:1,transpose(1./(.01.*Resilience(1:3,:))))

** **