Water Quality Modeling

Developing a GIS application to evaluate nitrogen and phosphorus land surface loads for the Nueces Bay

by

Zsuzsanna Iwanicka

GIS in Water Resources Fall 2002

Summary

This project attempted to develop a user friendly GIS application for total land surface loads calculation using Visual Basic for Application (VBA) in the ArcGIS environment. The mathematical model incorporated two steps: (i) determining the depth of the surface runoff in the function of precipitation and land use type and (ii) calculating selected forms of nitrogen and phosphorus loads released to the Nueces Bay from its the immediate watershed. The model returns the results in a display window. This paper includes the applied algorithms, the VBA program as well as describes the steps that need to be undertaken to utilize this application in ArcMap.

Objectives

In watersheds where the largest pollution loads are released through surface runoff, adequate water quality protection often requires a good understanding of how land use practices and its changes may affect the amount of released surface pollution loads. One of the challenges in the cooperation between the scientific community and policy makers is, that tools and models develop by scientists are not always easily applicable by other professional groups. Since in the policy making process decision makers needs to analyze possible options, there is a need for tools that could be utilized out side of the academia. The main purpose of this project is to learn to what extent the ESRI ArcGIS software may be applied to develop a user friendly tool that could be relevant for determining the impact of land use changes on non-point source loads releases through a defined watershed.

Methodology

The mathematical model of surface loads was taken from Ann M. Quenzer’s Master Thesis “GIS Assessment of the Total Loads and Water Quality in the Corpus Christi Bay System”. The model runs on the assumptions, that surface runoff as well as pollution concentrations are the function of land use. In the first step the model calculates the mean annual runoff depth as a function precipitation and land use type. In the second step the model calculates five constituents load Total Nitrogen, Total Kjeldahl Nitrogen, Nitrate and Nitrite, Total Phosphorus and Dissolved Phosphorus. Applied equations and constituent concentrations listed in Table 1 were taken from Ann Quenzer’s master thesis.

Table 1.Runoff depth equations and Event Mean Concentration Values for selected land use types

Event Mean Concentration Values (EMC) for the Modeled Parameters	Urban	Agricultural	Range land
Runoff depth[mm/year]	Q=0.24*P	Q=0.008312*exp(0.11415*P)	Q=0.008312*exp(0.11415*P)
Total Nitrogen [gN/m^3]	1.57	4.4	0.7
Total Kjeldahl Nitrogen [gN/m^3]	1.25	1.7	0.2
Nitrate and Nitrite [gN/m^3]	0.34	1.6	0.01
Total Phosphorus [gP/m^3]	0.35	1.3	-
Dissolved Phosphorus [gP/m^3]	0.23	-	< 0.01

P – annual mean precipitation mm/year

The first step in the model includes the calculation of depth of the surface runoff for each land use type using the equations presented in the first row of table 1. The land surface loads were calculated for each of the modeled constituents using the following equation:

Load [kg/d] = Q [mm/year] * Concentration [g/m³] * Area [m²] *10^-3m/mm * 10^-3kg/g * 1yr/365.2d

GIS Data preparation

The first step of data preparation included the definition of the study area through watershed delineation. Next, the defined study area was clipped to the land use polygons, then to the annual precipitation polygons. The results of these steps are presented in Figure 1. Land use and precipitation data was downloaded from the Texas Natural Resources Information Systems (TNRIS) website. The delineated watershed was obtained from Imane Mirni, a graduate student and researcher at the University of Texas Center for Research in Water Resources (CRWR).

Figure 1. Land use and precipitation layers of the study area

In the next step the land use and the precipitation data layers were intersected using the Geoprocessing Wizard. It allowed generating and store all the needed information for each particular land use polygon such as area, land use type and annual rainfall values in a single attribute table. The result of the intersection is presented in on Figure 2.

Figure 2. Intersected layer of the study area and its attribute table

Visual Basic Application

The main reason for developing this model in VBA was that this approach would allow conducting the analysis within one single application such as ArcMap. The spatial nature of the land use planning suggest that using graphical queries maybe a more convenient way to select and modify land use types than making the same modifications in tabular data format especially that is not connected directly to spatial display. The developed VBA application (commands included in App.A) calculates the total non-point source loads for each of the defined constituents. In order to analyze the effect of land use change on land surface loads, the land use codes need to be changed. These codes are included in the “LULC” field in the attribute table. The modification of the land use codes can be done only if the editing mode is active. (View/Tool bars/Editor/ Start Editing - as a task choose “Modify Feature”). The changes of land use codes can be done two ways:

• selecting the area (using the standard ArcMap “Select Features” button) and opening attribute table and changing the land use code in the LULC field.
• using a VBA tool developed especially for this purpose: like in the first option it requires to select an area using the “Select Features” button and than click on the designed button (in red circle on figure 3). It executes a set of commends (included in App.B) that calls out the “Attributes window” and allows to modify the land use type in that window without opening the whole attribute table.

Figure 3: Example of using the attribute window to modify land use codes

Once the landuse types are determined the new loads can be calculated by clicking on the button in the blue oval. It activates the VBA commends and generates the final result of the application is presented on figure 4 below. This window offers the option of copying the values to a separate file

Figure 4. The result...

Conclusions

The final product of this project indicates that ArcGIS provides several opportunities for water quality modeling. Beside simple models, like the one developed in this project, complex and powerful models ideas could be explored, for example incorporating the ArcHydro data and command structure. The second objective of this project: developing a user friendly application can be considered as partially successful. Although the model runs within one application, it still requires ArcGIS or at least ArcMap knowledge, which in the present might limit broader application.

Special Thanks…

To my husband, Jan Iwanicki, for his help in developing the VBA commands

To Imane Mirni, for sharing her GIS data

To Professor Maidment for being a wonderful teacher

Appendix A: The VBA commends to calculate non point surface loads:

Private Sub CommandButton2_Click()

'*********************************************

'Runoff and load calculation on Feature Dataset

'*********************************************

On Error GoTo Shame_end

'objects

Dim pWorkspaceFactory As IWorkspaceFactory

Dim pWorkspace As IWorkspace

Dim pFeatureWorkspace As IFeatureWorkspace

Dim pTable As ITable

Dim pRow As IRow

Dim pCursor As ICursor

Dim pFeatureDataset As IFeatureDataset

Dim pEnumDataset As IEnumDataset

'vars

Dim iLandUse As Long

Dim iArea As Long

Dim iRainFall As Long

Dim iRunoff As Double

Dim iP As Double

Dim iPO4 As Double

Dim iTKN As Double

Dim iNO2_3 As Double

Dim iN1 As Double

Dim Sum As Variant

Dim Txt As String

Dim R As Long

Dim I As Long

Dim Area As Double

Dim Q As Double

Dim QTot As Double

Dim P As Double 'Total Phosphorus partial value by landuse type[kg/d]

Dim PTot As Double 'Total Phosphorus for the modeled basin [kg/d]

Dim PO4 As Double 'Dissolved Phosphorus partial value by landuse type[kg/d]

Dim PO4Tot As Double 'Dissolved Phosphorus for the modeled basin [kg/d]

Dim TKN As Double 'Total Kjeldahl Nitrogen partial value by landuse type[kg/d]

Dim TKNTot As Double 'Total Kjeldahl Nitrogen [kg/d]for the modeled basin [kg/d]

Dim N As Double 'Total Nitrogen partial value by landuse type[kg/d]

Dim NTot As Double 'Total Nitrogen load for the modeled basin [kg/d]

Dim NO2_3 As Double 'Nitrate and Nitrite load partial value by landuse type [kg/d]

Dim NO2_3Tot As Double 'Nitrate and Nitrite for the modeled area [kg/d]

Set pWorkspaceFactory = New AccessWorkspaceFactory

Set pWorkspace = pWorkspaceFactory.OpenFromFile("D:\tp_data\lu_rain_c_p.mdb", 0)

Set pFeatureWorkspace = pWorkspace

Set pTable = pFeatureWorkspace.OpenTable("lu_rain_c")

'File indexes

iArea = pTable.FindField("Shape_Area")

iLandUse = pTable.FindField("LULC")

iRainFall = pTable.FindField("RANGE")

iRunoff = pTable.FindField("Runoff")

iP = pTable.FindField("TOT_P_MGP_")

iPO4 = pTable.FindField("PO4_MGP_L")

iTKN = pTable.FindField("TKN_MGN_L")

iNO2_3 = pTable.FindField("NO2_NO3_MG")

iN1 = pTable.FindField("TOT_N_MGN_")

'Number of records

Dim pQueryFilter As IQueryFilter

Set pQueryFilter = New QueryFilter

R = pTable.RowCount(pQueryFilter)

'Calculations loop

Set pCursor = pTable.Search(Nothing, False)

I = 0

Set pRow = pCursor.NextRow

UserForm1.MousePointer = fmMousePointerHourGlass

Do Until pRow Is Nothing

I = I + 1

'Update record counter

TextBox2.Value = I & " of " & R

DoEvents

Area = pRow.Value(iArea)

'Runoff & Loads

Select Case pRow.Value(iLandUse)

Case 11 To 17 'landuse type: urban

Q = 0.24 * (pRow.Value(iRainFall)) * 0.0254 * Area

P = Q * 0.35 / 365250 'P = 0.35[g/m3]

PO4 = Q * 0.23 / 365250 'PO4 = 0.23[g/m3]

TKN = Q * 1.25 / 365250 'TKN = 1.25[g/m3]

NO2_3 = Q * 0.34 / 365250 'NO2_3 = 0.34[g/m3]

N = Q * 1.57 / 365250 'N = 1.57 [g/m3]

Case 21, 23, 24 'landuse type: agricultural

Q = 0.008312 * Exp(0.011415 * pRow.Value(iRainFall) * 0.0254) * Area

P = Q * 1.3 / 365250

PO4 = Q * 0 / 365250 'PO4 no value reported

TKN = Q * 1.7 / 365250

NO2_3 = Q * 1.6 / 365250

N = Q * 4.4 / 365250

Case 31, 32, 33, 75, 76 'landuse type: rangeland

Q = 0.0053 * Exp(0.010993 * pRow.Value(iRainFall) * 0.0254) * Area

P = Q * 0 / 365250 'P no value reported

PO4 = Q * 0 / 365250 'PO4 < 0.01

TKN = Q * 0.2 / 365250

NO2_3 = Q * 0.01 / 365250

N = Q * 0.7 / 365250

Case 51, 52, 53, 54, 61, 62 'landuse type: water

Q = pRow.Value(iRainFall) * 0.0254 * Area

P = 0

PO4 = 0

TKN = 0

NO2_3 = 0

N = 0

End Select

'Store in database

pRow.Value(iRunoff) = Q

pRow.Value(iP) = P

pRow.Value(iPO4) = PO4

pRow.Value(iTKN) = TKN

pRow.Value(iNO2_3) = NO2_3

pRow.Value(iN1) = N

pRow.Store

Set pRow = pCursor.NextRow

'Keep total values

QTot = QTot + Q

PTot = PTot + P

PO4Tot = PO4Tot + P

TKNTot = TKNTot + TKN

NO2_3Tot = NO2_3Tot + NO2_3

NTot = NTot + N

Loop

'Output total data to form

UserForm1.MousePointer = fmMousePointerDefault

Txt = ""

Txt = Txt & " Q = " & Format(QTot, "##,##0.00") & " m3/yr" & Chr(13) & Chr(10)

Txt = Txt & " P = " & Format(PTot, "##,##0.00") & " kgP/d" & Chr(13) & Chr(10)

Txt = Txt & " P04 = " & Format(PO4Tot, "##,##0.00") & " kgP/d" & Chr(13) & Chr(10)

Txt = Txt & " TKN = " & Format(TKNTot, "##,##0.00") & " kgN/d" & Chr(13) & Chr(10)

Txt = Txt & " NO2/3 = " & Format(NO2_3Tot, "##,##0.00") & " kgN/d" & Chr(13) & Chr(10)

Txt = Txt & " N = " & Format(NTot, "##,##0.00") & " kgN/d" & Chr(13) & Chr(10)

TextBox1.Value = Txt

Exit Sub

Shame_end:

TextBox1.Value = "!!! Calculation or Data Error !!!"

End Sub

Private Sub CommandButton3_Click()

'*****************************************

'Simple copy to clipboard

'*****************************************

TextBox1.SetFocus

TextBox1.SelStart = 0

TextBox1.SelLength = TextBox1.TextLength

SendKeys ("^c")

End Sub

Appendix B.VBA Commends to change the land use code through the attribute window

Public Sub ToggleAttributeWindow()

'Shows getting the Editor's reference by class ID

Dim pEditor As IEditor

Dim pID As New UID

pID = "esriCore.Editor"

Set pEditor = Application.FindExtensionByCLSID(pID)

Dim pExtension As IExtension

Dim pExtensionManager As IExtensionManager

Dim pAttributeWindow As IAttributeWindow

Dim iExtensionCount As Integer

'Ensure an active edit session

If Not pEditor.EditState = esriStateEditing Then

MsgBox "Turn to 'EDIT MODE' first."

Exit Sub

End If

'Loop through all of the extensions to find the attribute window

Set pExtensionManager = pEditor 'QI

For iExtensionCount = 0 To pExtensionManager.ExtensionCount - 1

Set pExtension = pExtensionManager.Extension(iExtensionCount)

If pExtension.Name = "AttributeWindow" Then Exit For

Next iExtensionCount

'Open the window if it is closed; close it, if it is open

If Not TypeOf pExtension Is IAttributeWindow Then Exit Sub

Set pAttributeWindow = pExtension 'QI

If Not pAttributeWindow.Visible Then

pAttributeWindow.Visible = True

Else

pAttributeWindow.Visible = False

End If

End Sub