CE 470: Stream Riparian Systems and Fluvial Morphology

Name: Cassio Rampinelli

Date: September 26, 2019

Assignment: Homework 2

Note: This script was done in R programming language

______________________________________________________________________________

1.PURPOSE

This assignment presents a one-dimensional hydraulic model to compute water surface profile for irregular cross-sections. The code was developed in R programming language in such way that can be easily adjusted to any number and geometry of cross-sections. In order to evaluate the quality of the developed model, the results were compared to an HEC-RAS model output by employing equivalent settings. The model focus on backwater curves that normally occur upstream of control structures such as dams and weirs. Such models are useful to assess the water level upstream of hydraulic structures and consequently areas that can be eventually flooded after the implementation of such structures.

The methodology, assumptions and basic equations are presented in the next section followed by the results and and discussion.

2.METHODS

2.1 Theoretical Background

A basic assumption in calculating water-surface profiles is that the head loss between upstream and downstream sections can be estimated using Manning’s equation to estimate the slope of the energy grade line ($S_f$), given by:

\[ S_f=\Bigg(\frac{n\cdot Q}{1.49\cdot A\cdot R^{2/3}}\Bigg)^2 \] In which

$S_f$= Slope of the energy gradeline

$Q$= Flow ($cfs$)

$A$= Flow area ($ft^2$)

$R$ = Hydraulic Radius ($ft$)

In this way it is assumed that flow conditions change gradually.Water surface profiles are computed from one cross section to the next by solving the Energy Equation with an iterative procedure. The procedure can use: the Direct-integration method, or the Direct-step method, or the Standard-step method. In this work the Standard-step method was emplyoed. Figure 1 presents the terms related to the Energy Equation that is written as follows:

\[ E_1=E_2 \]

\[ WS_1+\frac{V^2_1}{2g}+h_L=WS_2+\frac{V^2_2}{2g} \]

Figure 2.1. Representation of the terms related to the Energy Equation

The Energy Head Loss ($h_L$) between two cross sections is comprised of friction losses and contraction or expansion losses. The equation for the energy head loss is as follows.

\[ h_L=\bar{S_f}L+C\Bigg|\alpha_2\frac{V_2^2}{2g}-\alpha_1\frac{V_1^2}{2g}\Bigg| \]

In which (in English units)

$h_f$= head loss ($ft$)

$\bar{S_f}$= Average gradeline slope

$C$= Expansion (0.3) or contraction (0.1) loss coefficient

$\alpha_1, \alpha_2$= velocity weighting coefficients. (Assumed to be 1 in this work)

$V, V_2$= Average velocities ($ft/s$) $nbsp

g=gravitational acceleration (32.2 $ft/s^2$)

$R$ = Hydraulic Radius ($ft$)

2.2 Standard Step Method Routine

The Standard Step Method routine can be summarized as follows:

Make a gess for the water surface elevation $WS_2$ for the upstream section (an initial gess can be the same $WS_1$);
Compute ($E_2$-$E_1$) that is the head loss based on the Energy Equation;
Compute the head loss $h_L$ based on the Energy Head Loss Equation;
Compare (b) and (c). If no match, estimate $WS_2$ using the lead loss Equation. Then, compute a new $WS_2$ using the average between the new $WS_2$ and the previous one;
Repeat the procedure until the difference between the previous and the new $WS_2$ is acceptable. (In this simulation we assumed a difference less than 0.01 reasonable ).

2.3 Data

The boundary conditions assumed for the simulation was a flow rate of 100 cfs and a downstream depth of 5’. The used cross-sections are presented in the tables below. For simplicity Manning values were attributed to the left bank, channel and right bank of each cross-section

Figure 2.2 shows a plot of all cross-sections together. The cross-sections are equidistant of 2000 fts.

nbsp;

Table 2.1. Cross Sections

The Cross Sections Plots are presented bellow. The dots indicate the partition between the channel and banks.

nbsp;

x1<-c(0,3,10,18,23,31,38,41)
y1<-c(7,4,4,0,0,4,4,7)

x2<-c(1,5,15,20,28,30,35)
y2<-c(9,4,4,2,2,4,9)

x3<-c(0,4,12,17,25,30,38,38.1)
y3<-c(11,8,8,4,4,9,9,11)

plot(x1,y1,type="l",main="Cross-Sections", xlab="Station(ft)",ylab="Elevation(f)",xlim=c(0,40),ylim=c(0,12),lwd=3)
lines(x2,y2,lty=2,lwd=3,col="red")
lines(x3,y3,lty=3,lwd=2,col="darkgreen")
grid()
legend("top",lty=c(1,2,3),lwd=c(3,3,2),bty="n",cex=1.5,col=c("black","red","darkgreen"),
       legend=c("Sec.1","Sec.2","Sec.3"))

2.4 Code

The implemented code is presented as follows with some comments. The code basically has 2 auxiliary functions (wettedArea.Perim and deltaEcalc) and a main function (findWsd). The first function computes the wetted area and perimeter. The second one computes the head loss and the last one is the main function that computes the water surface elevation based on the standard step method routine.

Figure 2.2 shows a summary of the functions developed in the code. All inputs and outputs of each function are specified in code line commentaries. At the end of the code the final water surface levels are shown followed by the geometric parameters for each section.

Figure 2.2. Functions developed in the code

nbsp;

#Clear memory
rm(list=ls(all=TRUE)) 


#Cross Sections Coordinates

x1<-c(0,3,10,18,23,31,38,41)
y1<-c(7,4,4,0,0,4,4,7)

x2<-c(1,5,15,20,28,30,35)
y2<-c(9,4,4,2,2,4,9)

x3<-c(0,4,12,17,25,30,38,38.1)
y3<-c(11,8,8,4,4,9,9,11)



#################
#Function to compute wetted area and wetted perimeter for anny geometry
###################

wettedArea.Perim<-function(x,y,y0,xleft,yleft,xright,yright)
{
  
  ##INPUTS
  #x= vector with x coordinates
  #y= vector with y coordinates
  #y0=water surface elevation
  #xleft=x left bank
  #yleft=y left bank
  #xright=x right bank
  #yright=y right bank
  
  ##OUTPUT: a List with 8 elements
  ##[[1]]Left Bank Wetted Area
  ##[[2]]Channel Wetted Area
  ##[[3]]Right Bank Wetted Area
  ##[[4]]Total Wetted Area
  ##[[5]]Left Bank Wetted Perimeter
  ##[[6]]Channel Wetted Perimeter
  ##[[7]]Right Bank Wetted Perimeter
  ##[[8]]Total Wetted Perimeter
  
  #########
  
  x.left=xleft
  y.left=yleft
  x.right=xright
  y.right=yright
  
  
  
  if(y0<=y.right && y0<=y.left)
  {
    
    
    x.left.side=x[1:ceiling(length(x)/2)]
    x.channel.left.side=x.left.side[max(which(x.left.side==x.left)):length(x.left.side)]
    
    y.left.side=y[1:ceiling(length(y)/2)]
    y.channel.left.side=y.left.side[(which(x.left.side>=x.left))]
    
    x.right.side=x[(ceiling(length(x)/2)+1):length(y)]
    x.channel.right.side=x.right.side[(which(x.right.side<=x.right))]
    
    y.right.side=y[(ceiling(length(x)/2)+1):length(y)]
    y.channel.right.side=y.right.side[(which(x.right.side<=x.right))]
    
    
    x0left=x.channel.left.side[length(x.channel.left.side)]-(y0- y.channel.left.side[length(y.channel.left.side)])*(x.channel.left.side[length(x.channel.left.side)]-x.left)/(y.left-y.channel.left.side[length(y.channel.left.side)])
    
    
    x0right=x.channel.right.side[min(which(x.channel.right.side<x.right))]+(y0-y.channel.right.side[1])*(x.right-x.channel.right.side[min(which(x.channel.right.side<x.right))])/(y.right-y.channel.right.side[1])
    
    
    #Computing channel area
    xx=yy=0
    xx=c(x0left,x.channel.left.side[which(x.channel.left.side>x0left)],x.channel.right.side[which(x.channel.right.side<x0right)],x0right)
    xx[2:length(xx)]=xx[length(xx):2]
    
    yy=c(y0,y.channel.left.side[which(x.channel.left.side>x0left)],y.channel.right.side[which(x.channel.right.side<x0right)],y0)
    yy[2:length(yy)]=yy[length(yy):2]
    
    a.channel=sum((sum(xx[1:length(xx)-1]*yy[2:length(xx)])+xx[length(xx)]*yy[1])-sum(xx[2:length(xx)]*yy[1:length(xx)-1])-xx[1]*yy[length(xx)])
    a.channel=abs(a.channel)/2
    
    
    #Computing perimeter
    
    xxp=yyp=0
    xxp=c(x0left,x.channel.left.side[which(x.channel.left.side>x0left)],x.channel.right.side[which(x.channel.right.side<x0right)],x0right)
    yyp=c(y0,y.channel.left.side[which(x.channel.left.side>x0left)],y.channel.right.side[which(x.channel.right.side<x0right)],y0)
    
    xsquared=(xxp[1:length(xxp)-1]-xxp[2:length(xxp)])^2
    ysquared=(yyp[1:length(yyp)-1]-yyp[2:length(yyp)])^2
    
    dists=sqrt(xsquared+ysquared)
    perim.channel=sum(dists)
    perim.channel
    
    A.right=A.left=perim.left=perim.right=0
    
    
  }
  
  
  else{
    
    
    
    if(y0>y.left && y0>y.right)
    {
      
      
      #Computing left area
      
      x.left.side=x[1:ceiling(length(x)/2)]
      x.channel.left.side=x.left.side[max(which(x.left.side==x.left)):length(x.left.side)]
      x.bank.left.side=x.left.side[1:max(which(x.left.side==x.left))]
      
      y.left.side=y[1:ceiling(length(y)/2)]
      y.channel.left.side=y.left.side[max(which(x.left.side==x.left)):length(x.left.side)]
      y.bank.left.side=y.left.side[1:max(which(x.left.side==x.left))]
      
      x0left= x.bank.left.side[min(which(y.bank.left.side<y0))]-(y0-y.bank.left.side[min(which(y.bank.left.side<y0))])*(x[min(which(y.bank.left.side<y0))]-x.bank.left.side[match(min(y.bank.left.side[which(y.bank.left.side>y0)]),y.bank.left.side)])/(y.bank.left.side[match(min(y.bank.left.side[which(y.bank.left.side>y0)]),y.bank.left.side)]-y.bank.left.side[min(which(y.bank.left.side<y0))])
      
      xtemp.left=x.bank.left.side[which(x.bank.left.side>x0left)]
      xtemp.left=c(x0left,xtemp.left,xtemp.left[length(xtemp.left)])
      
      ytemp.left=y.bank.left.side[which(x.bank.left.side>x0left)]
      ytemp.left=c(y0,ytemp.left,y0)
      
      #Re-ording the x coordinates
      xx=0
      xx=xtemp.left
      xx[2:length(xx)]=xx[length(xx):2]
      
      #Re-ording the y coordinates
      yy=0
      yy=ytemp.left
      yy[2:length(yy)]=yy[length(yy):2]
      
      #Computing the area below the water surface level
      A.left=sum((sum(xx[1:length(xx)-1]*yy[2:length(xx)])+xx[length(xx)]*yy[1])-sum(xx[2:length(xx)]*yy[1:length(xx)-1])-xx[1]*yy[length(xx)])
      A.left=abs(A.left)/2
      
      #Computing perimeter
      xxp=xtemp.left[1:(length(xtemp.left)-1)]
      yyp=ytemp.left[1:(length(ytemp.left)-1)]
      
      xsquared=(xxp[1:length(xxp)-1]-xxp[2:length(xxp)])^2
      ysquared=(yyp[1:length(yyp)-1]-yyp[2:length(yyp)])^2
      
      dists=sqrt(xsquared+ysquared)
      perim.left=sum(dists)
      
      
      
      #Computing right area
      x.right.side=x[(ceiling(length(x)/2)+1):length(x)]
      
      x.channel.right.side=x.right.side[1:max(which(x.right.side==x.right))]
      
      x.bank.right.side=x.right.side[max(which(x.right.side==x.right)):length(x.right.side)]
      
      y.right.side=y[length(y):(ceiling(length(y)/2)+1)]
      y.right.side=y.right.side[length(y.right.side):1]
      
      y.channel.right.side=y.right.side[1:max(which(x.right.side==x.right))]
      
      y.bank.right.side=y.right.side[length(y.right.side):max(which(x.right.side==x.right))]
      y.bank.right.side=y.bank.right.side[length(y.bank.right.side):1]
      
      x0right= x.bank.right.side[max(which(y.bank.right.side<=y0))]+(y0-y.bank.right.side[max(which(y.bank.right.side<=y0))])*(x.bank.right.side[min(which(y.bank.right.side>y0))]-x.bank.right.side[max(which(y.bank.right.side<=y0))])/(y.bank.right.side[min(which(y.bank.right.side>y0))]-y.bank.right.side[max(which(y.bank.right.side<=y0))])
      
      xtemp.rigth=x.bank.right.side[which(x.bank.right.side<x0right)]
      
      xtemp.rigth=c(xtemp.rigth[1],xtemp.rigth,x0right)
      
      
      ytemp.right=y.bank.right.side[which(x.bank.right.side<x0right)]
      ytemp.right=c(y0,ytemp.right,y0)
      
      #Re-ording the x coordinates
      xx=0
      xx=xtemp.rigth
      xx[2:length(xx)]=xx[length(xx):2]
      
      #Re-ording the y coordinates
      yy=0
      yy=ytemp.right
      yy[2:length(yy)]=yy[length(yy):2]
      
      #Computing the area below the water surface level
      A.right=sum((sum(xx[1:length(xx)-1]*yy[2:length(xx)])+xx[length(xx)]*yy[1])-sum(xx[2:length(xx)]*yy[1:length(xx)-1])-xx[1]*yy[length(xx)])
      A.right=abs(A.right)/2
      
      #Computing perimeter
      
      xxp=yyp=0
      xxp=xtemp.rigth[2:(length(xtemp.rigth))]
      yyp=ytemp.right[2:(length(ytemp.right))]
      
      xsquared=(xxp[1:length(xxp)-1]-xxp[2:length(xxp)])^2
      ysquared=(yyp[1:length(yyp)-1]-yyp[2:length(yyp)])^2
      
      dists=sqrt(xsquared+ysquared)
      perim.right=sum(dists)
      perim.right
      
      #Computing channel area
      xx=yy=0
      xx=c(x.channel.left.side[1],x.channel.left.side,x.channel.right.side,x.channel.right.side[length(x.channel.right.side)])
      xx[2:length(xx)]=xx[length(xx):2]
      
      yy=c(y0,y.channel.left.side,y.channel.right.side,y0)
      yy[2:length(yy)]=yy[length(yy):2]
      
      a.channel=sum((sum(xx[1:length(xx)-1]*yy[2:length(xx)])+xx[length(xx)]*yy[1])-sum(xx[2:length(xx)]*yy[1:length(xx)-1])-xx[1]*yy[length(xx)])
      a.channel=abs(a.channel)/2
      
      
      #Computing perimeter
      
      xxp=yyp=0
      xxp=c(x.channel.left.side,x.channel.right.side)
      yyp=c(y.channel.left.side,y.channel.right.side)
      
      xsquared=(xxp[1:length(xxp)-1]-xxp[2:length(xxp)])^2
      ysquared=(yyp[1:length(yyp)-1]-yyp[2:length(yyp)])^2
      
      dists=sqrt(xsquared+ysquared)
      perim.channel=sum(dists)
      
      
      
    }
    
    
    else
    {
      
      
      if(y0>y.left && y0<=y.right)
      {
        
        #Computing left area
        
        x.left.side=x[1:ceiling(length(x)/2)]
        x.channel.left.side=x.left.side[max(which(x.left.side==x.left)):length(x.left.side)]
        x.bank.left.side=x.left.side[1:max(which(x.left.side==x.left))]
        
        y.left.side=y[1:ceiling(length(y)/2)]
        y.channel.left.side=y.left.side[max(which(x.left.side==x.left)):length(x.left.side)]
        y.bank.left.side=y.left.side[1:max(which(x.left.side==x.left))]
        
        x0left= x.bank.left.side[min(which(y.bank.left.side<y0))]-(y0-y.bank.left.side[min(which(y.bank.left.side<y0))])*(x[min(which(y.bank.left.side<y0))]-x.bank.left.side[match(min(y.bank.left.side[which(y.bank.left.side>y0)]),y.bank.left.side)])/(y.bank.left.side[match(min(y.bank.left.side[which(y.bank.left.side>y0)]),y.bank.left.side)]-y.bank.left.side[min(which(y.bank.left.side<y0))])
        
        xtemp.left=x.bank.left.side[which(x.bank.left.side>x0left)]
        xtemp.left=c(x0left,xtemp.left,xtemp.left[length(xtemp.left)])
        
        ytemp.left=y.bank.left.side[which(x.bank.left.side>x0left)]
        ytemp.left=c(y0,ytemp.left,y0)
        
        #Re-ording the x coordinates
        xx=0
        xx=xtemp.left
        xx[2:length(xx)]=xx[length(xx):2]
        
        #Re-ording the y coordinates
        yy=0
        yy=ytemp.left
        yy[2:length(yy)]=yy[length(yy):2]
        
        #Computing the area below the water surface level
        A.left=sum((sum(xx[1:length(xx)-1]*yy[2:length(xx)])+xx[length(xx)]*yy[1])-sum(xx[2:length(xx)]*yy[1:length(xx)-1])-xx[1]*yy[length(xx)])
        A.left=abs(A.left)/2
        
        #Computing perimeter
        xxp=xtemp.left[1:(length(xtemp.left)-1)]
        yyp=ytemp.left[1:(length(ytemp.left)-1)]
        
        xsquared=(xxp[1:length(xxp)-1]-xxp[2:length(xxp)])^2
        ysquared=(yyp[1:length(yyp)-1]-yyp[2:length(yyp)])^2
        
        dists=sqrt(xsquared+ysquared)
        perim.left=sum(dists)
        
        A.right=perim.right=0
        
        #Computing channel area
        xx=yy=0
        xx=c(x.channel.left.side[1],x.channel.left.side,x.channel.right.side)
        xx[2:length(xx)]=xx[length(xx):2]
        
        yy=c(y0,y.channel.left.side,y.channel.right.side)
        yy[2:length(yy)]=yy[length(yy):2]
        
        a.channel=sum((sum(xx[1:length(xx)-1]*yy[2:length(xx)])+xx[length(xx)]*yy[1])-sum(xx[2:length(xx)]*yy[1:length(xx)-1])-xx[1]*yy[length(xx)])
        a.channel=abs(a.channel)/2
        
        
        #Computing perimeter
        
        xxp=yyp=0
        xxp=c(x.channel.left.side,x.channel.right.side)
        yyp=c(y.channel.left.side,y.channel.right.side)
        
        xsquared=(xxp[1:length(xxp)-1]-xxp[2:length(xxp)])^2
        ysquared=(yyp[1:length(yyp)-1]-yyp[2:length(yyp)])^2
        
        dists=sqrt(xsquared+ysquared)
        perim.channel=sum(dists)
        perim.channel
        
        
        
        
        
      }
      
      else
      {
        
        if(y0<y.left && y0>y.right)
        {
          
          
          #Computing right bank area
          x.right.side=x[(ceiling(length(x)/2)+1):length(x)]
          
          x.channel.right.side=x.right.side[1:max(which(x.right.side==x.right))]
          
          x.bank.right.side=x.right.side[max(which(x.right.side==x.right)):length(x.right.side)]
          
          y.right.side=y[length(y):(ceiling(length(y)/2)+1)]
          y.right.side=y.right.side[length(y.right.side):1]
          
          y.channel.right.side=y.right.side[1:max(which(x.right.side==x.right))]
          
          y.bank.right.side=y.right.side[length(y.right.side):max(which(x.right.side==x.right))]
          y.bank.right.side=y.bank.right.side[length(y.bank.right.side):1]
          
          x0right= x.bank.right.side[max(which(y.bank.right.side<=y0))]+(y0-y.bank.right.side[max(which(y.bank.right.side<=y0))])*(x.bank.right.side[min(which(y.bank.right.side>y0))]-x.bank.right.side[max(which(y.bank.right.side<=y0))])/(y.bank.right.side[min(which(y.bank.right.side>y0))]-y.bank.right.side[max(which(y.bank.right.side<=y0))])
          
          xtemp.rigth=x.bank.right.side[which(x.bank.right.side<x0right)]
          
          xtemp.rigth=c(xtemp.rigth[1],xtemp.rigth,x0right)
          
          
          ytemp.right=y.bank.right.side[which(x.bank.right.side<x0right)]
          ytemp.right=c(y0,ytemp.right,y0)
          
          #Re-ording the x coordinates
          xx=0
          xx=xtemp.rigth
          xx[2:length(xx)]=xx[length(xx):2]
          
          #Re-ording the y coordinates
          yy=0
          yy=ytemp.right
          yy[2:length(yy)]=yy[length(yy):2]
          
          #Computing the area below the water surface level
          A.right=sum((sum(xx[1:length(xx)-1]*yy[2:length(xx)])+xx[length(xx)]*yy[1])-sum(xx[2:length(xx)]*yy[1:length(xx)-1])-xx[1]*yy[length(xx)])
          A.right=abs(A.right)/2
          
          #Computing perimeter
          
          xxp=yyp=0
          xxp=xtemp.rigth[2:(length(xtemp.rigth))]
          yyp=ytemp.right[2:(length(ytemp.right))]
          
          xsquared=(xxp[1:length(xxp)-1]-xxp[2:length(xxp)])^2
          ysquared=(yyp[1:length(yyp)-1]-yyp[2:length(yyp)])^2
          
          dists=sqrt(xsquared+ysquared)
          perim.right=sum(dists)
          
          A.left=perim.left=0
          
          
          #Computing channel area
          xx=yy=0
          xx=c(x.channel.left.side,x.channel.right.side,x.channel.right.side[length(x.channel.right.side)])
          xx[2:length(xx)]=xx[length(xx):2]
          
          yy=c(y.channel.left.side,y.channel.right.side,y0)
          yy[2:length(yy)]=yy[length(yy):2]
          
          a.channel=sum((sum(xx[1:length(xx)-1]*yy[2:length(xx)])+xx[length(xx)]*yy[1])-sum(xx[2:length(xx)]*yy[1:length(xx)-1])-xx[1]*yy[length(xx)])
          a.channel=abs(a.channel)/2
          
          
          #Computing perimeter
          
          xxp=yyp=0
          xxp=c(x.channel.left.side,x.channel.right.side)
          yyp=c(y.channel.left.side,y.channel.right.side)
          
          xsquared=(xxp[1:length(xxp)-1]-xxp[2:length(xxp)])^2
          ysquared=(yyp[1:length(yyp)-1]-yyp[2:length(yyp)])^2
          
          dists=sqrt(xsquared+ysquared)
          perim.channel=sum(dists)
          perim.channel
          
          
          
          
        }
        
      }
      
    }
    
  }
  
  
  a.total=A.left+a.channel+A.right
  perim.total=perim.left+perim.channel+perim.right
  
  result=list(A.left,a.channel,A.right,a.total,perim.left,perim.channel,perim.right,perim.total)
  
  return(result) 
  
}


#################
#Function to compute the head loss
###################


deltaEcalc<-function(sfd,sfu,vd,vu,dx){
  
  Sf1=sfd
  Sf2=sfu
  
  deltax=dx
  v1=vd
  v2=vu
  
  Sfbar=mean(c(Sf1,Sf2))
  
  veld=vd
  velu=vu
  
  if(veld>velu){C=0.1}
  else C=0.3
  
  #C=0.3
  
  
  deltaEcalc=Sfbar*deltax+C*abs((velu^2)/(2*32.2)-(veld^2)/(2*32.2))
  
  return(deltaEcalc)
}




###############################
#Standard step method function
###############################


findWsd<-function(ws1,v1,E1,Sf1,xup,yup,yd,posD,posU,nL,nC,nR,xleft,yleft,xright,yright){
                             
  
  ##INPUTS
  #ws1= water surface elevation from downstream section
  #v1= velocity from downstream section
  #E1=Energy from downstream section
  #Sf1= Slope of the Energy grade line at downstream section
  #xup=x coordinates from the upstream section
  #yup=y coordinates from the upstream section
  #yd= y coordinates from the downstream section
  #posD=Station of the downstream section
  #posU=Station of the upstream section
  #nL=Left bank Manning
  #nC=Channel Manning
  #nR-Right bank Manning
  #xleft=x left bank
  #yleft=y left bank
  #xright=x right bank
  #yright=y right bank
  
  ##OUTPUT: Water surface elevation from the downstream section
  
  ws1=ws1
  v1=v1
  E1=E1
  Sf1=Sf1
  xup=xup
  yup=yup
  yd=yd
  posD=posD
  posU=posU
nL=nL
nC=nC
nR=nR
xleft=xleft
yleft=yleft
xright=xright
yright=yright
  

  counter=1
  wsu=(ws1-min(yd))+min(yup)
  
  
  while(counter==1){
    
    
    #Section 2
    #y02=min(y2)+(y0-min(y1))
    ws2=wsu
    
    restult=wettedArea.Perim(x=xup,y=yup,y0=ws2,xleft=xleft,yleft=yleft,xright=xright,yright=yright)
    A2=restult[[4]]
    v2=Q/A2
    P2=restult[[8]]
    Rh2=A2/P2
    E2=min(yup)+(ws2-min(yup))+v2^2/(2*32.2)
    
    PL=restult[[5]]
    PC=restult[[6]]
    PR=restult[[7]]
    
    n=((PL*(nL)^(1.5)+PC*(nL)^(1.5)+PR*(nR)^(1.5))/P2)^(2/3)
    
    Sf2=((n*Q)/(1.49*A2*(Rh2^(2/3))))^2
    
    #Checking
    deltaE=E2-E1
    
    deltaE.calc=deltaEcalc(Sf1,Sf2,v1,v2,(stations[posU]-stations[posD]))
    
    
    ws2.calc=ws1+(v1^2)/(2*32.2)+deltaE.calc-((v2^2)/(2*32.2))
    
    ws2.updated=(ws2.calc+ws2)/2
    
    delta=abs(ws2.updated-ws2)
    
    if(delta<0.001){
      
      counter=0
      return(ws2.updated)
    }
    
    else
      counter=1
    wsu=ws2.updated
    
    
  }
  
  
  
}


#####################################
#####################################
##RUNNING THE CODE##################
####################################


#Saving x an y matrix with the coordinates
xMatrix=qpcR:::cbind.na(x1,x2,x3)
yMatrix=qpcR:::cbind.na(y1,y2,y3)


#Stations of the sections
stations=c(0,2000,4000)

#Discharge (cfs)
Q=100

#Left Banks
xsleft<-c(x1[3],x2[3],x3[3])
ysleft<-c(y1[3],y2[3],y3[3])

#Right Banks
xsright<-c(x1[6],x2[6],x3[6])
ysright<-c(y1[6],y2[6],y3[6])

#Manning left bank
Nleft<-c(0.04,0.045,0.04)
#Manning channel
Nchannel<-c(0.025,0.024,0.025)
#Manning right bank
Nright<-c(0.035,0.024,0.035)



#Water surface elevation at downstream section(ft)
wsd=ws1=5
#Vector with the water surface elevation of the sections
WSE=ws1


#Initiating Variables
A=0
v=0
P=0
Rh=0
E=0
Sf=0
nsec=3


#downstream x coordinates without NA
x1<-xMatrix[,1]
x1<-x1[!is.na(x1)]

#downstream y coordinates without NA
y1<-yMatrix[,1]
y1<-y1[!is.na(y1)]


#Section downstream

restult=wettedArea.Perim(x=x1,y=y1,y0=ws1,xleft=xsleft[1],yleft=ysleft[1],xright=xsright[1],yright=ysright)

A[1]=restult[[4]]
v[1]=Q/A[1]
P[1]=restult[[8]]
Rh[1]=A[1]/P[1]
E[1]=wsd+v[1]^2/(2*32.2)

PL=restult[[5]]
PC=restult[[6]]
PR=restult[[7]]
n=((PL*(Nleft[1])^(1.5)+PC*(Nchannel[1])^(1.5)+PR*(Nright[1])^(1.5))/P[1])^(2/3)
Sf[1]=((n*Q)/(1.49*A[1]*(Rh[1]^(2/3))))^2


######################################
#Loop to apply the standard step method along the cross sections
######################################

for(i in 2:nsec){
  
  xd=yup=0
  
  #downstream x coordinates without NA
  xd<-xMatrix[,i-1]
  xd<-xd[!is.na(xd)]
  
  #downstream y coordinates without NA
  yd<-yMatrix[,i-1]
  yd<-yd[!is.na(yd)]
  
  #upwnstream x coordinates without NA
  xup<-xMatrix[,i]
  xup<-xup[!is.na(xup)]
  
  #upwnstream y coordinates without NA
  yup<-yMatrix[,i]
  yup<-yup[!is.na(yup)]
  
  
  wsup=findWsd(ws1=WSE[i-1],v1=v[i-1],E1=E[i-1],Sf1=Sf[i-1],xup=xup,yup=yup,yd=yd,posD=(i-1),posU=i,nL=Nleft[i],nC=Nchannel[i],nR=Nright[i],xleft=xsleft[i],yleft=ysleft[i],xright=xsright[i],yright=ysright[i])
  
  
  
  restult=wettedArea.Perim(x=xup,y=yup,y0=wsup,xleft=xsleft[i],yleft=ysleft[i],xright=xsright[i],yright=ysright[i])
  
  A[i]=restult[[4]]
  v[i]=Q/A[i]
  P[i]=restult[[8]]
  Rh[i]=A[i]/P[i]
  E[i]=wsd+v[i]^2/(2*32.2)
  
  PL=restult[[5]]
  PC=restult[[6]]
  PR=restult[[7]]
  n=((PL*(Nleft[i])^(1.5)+PC*(Nchannel[i])^(1.5)+PR*(Nright[i])^(1.5))/P[i])^(2/3)
  Sf[i]=((n*Q)/(1.49*A[i]*(Rh[i]^(2/3))))^2
  
  
  
  WSE[i]=wsup
  
  
  
  
}

#Return the water surface elevation for each section
WSE

## [1] 5.000000 5.767607 7.485919

#Return the wetted area for each cross section
A

## [1] 88.00000 70.00216 41.55794

#Return the velocity for each cross section
v

## [1] 1.136364 1.428527 2.406279

#Return the wetted perimeter for each cross section
P[i]=restult[[8]]
#Return the hydraulic radius for each cross section
Rh

## [1] 2.215678 2.259810 2.245158

#Return the Energy gradeline for each cross section
E

## [1] 5.020052 5.031688 5.089910

#Return the Energy slope
E

## [1] 5.020052 5.031688 5.089910

Sf

## [1] 0.0001888902 0.0003398438 0.0005544601

3.DISCUSSION

The plot below presents a comparison between the water surface profile from the code and from a run on HEC-RAS by applying the same settings. Although the computed water profile seems to be reasonable, it can be seen that it tends to overstemitate HEC-RAS model water profile.

plot(stations,c(min(y1),min(y2),min(y3)),type="l",main="Water Profile", xlab="Station(ft)",ylab="Elevation(f)",xlim=c(0,4000),ylim=c(0,7.5),lwd=3)

hec=c(5,5.31,6.81)
code=c(5,5.77,7.49)

lines(stations,hec,lty=2,lwd=3,col="red")
lines(stations,code,lty=1,lwd=3,col="blue")
legend("bottomright",lty=c(1,2,1),lwd=c(3,3,2),bty="n",cex=1.5,col=c("black","red","blue"),
       legend=c("Bed","W.S. HEC-RAS","W.S. Code"))
grid()

The Table below shows the difference between both outputs considering the velocity and the energy slope, besides the water profile. For the water surface elevation, on section 2 the proposed model overestimates in 8.7 % the water level predicted by HEC-RAS. On section 3 the difference increases and reaches 15.4%. The average error was 12% along the profile.

Considering the velocity, the discrepancies between the models were more significant with a general trend of the proposed model underestimate the velocity. The average error was -28.5% with regards the HEC-RAS computed velocities. The Energy Slope generated a variety of errors. For section 3 the energy slope computed by the code was more than 70% less than the one indicated by HEC-RAS output. However, for section 2 the difference between the models were only 3%. This parameter is on a small scale and consequently comparison between the results may be more sensitive to variations.

Finally, although some discrepancies were observed between the HEC-RAS model and the implemented one, the computed water surface profile can be considered reasonable. Part of the observed differences may be related to the omission of the velocity weighting coefficients ($\alpha1$ and $alpha2$) in the implemented code. Since the velocities between the sections had some considerable difference, and the HEC-RAS model takes these coefficients into consideration, it can increase the energy head loss what may lead to lower water surface elevation. However, more simulations and tests must be conducted to validate the implemented model.

$*Table 3.1. Comparison of the outputs*$

Table 3.1. Comparison of the outputs