diff --git a/.DS_Store b/.DS_Store
index 9c6f8ac..3ebc81f 100644
Binary files a/.DS_Store and b/.DS_Store differ
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..5e436a7
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,3 @@
+results_facility_level-*_HUC6.xlsx
+results_all_level-*_HUC6.pkl
+Data Dictionary.docx
\ No newline at end of file
diff --git a/README.md b/README.md
index f9ad818..c82ae80 100644
--- a/README.md
+++ b/README.md
@@ -17,6 +17,7 @@ One of the key challenges in the adoption of green infrastructure is the inheren
 2. **Data Input**: All inputs needed to run the code are found in the `inputs` folder. Input assumptions for the Green or Gray treatment methods can be updated in the `model_inputs.xlsx` spreadsheet.
 3. **Execution**: Run the script using `python green_wastewater_treatment_point-source_data-reduction.py`.
 4. **Results**: The script will output comparative results in terms of costs and emissions for each of the technologies which can be found in the `results` folder.
+5. **Post Processing**: After the results are generated, the `csv_creation_for_images.py` file can be used to generate csv files that summarize the results and can be used for figure generation.
 
 ## Contributions and Feedback
 We welcome contributions and feedback to improve this analysis. Please feel free to fork the repository, make changes, and submit a pull request. For feedback and suggestions, please open an issue in the repository.
diff --git a/csv_creation_for_images.py b/csv_creation_for_images.py
new file mode 100644
index 0000000..01968bf
--- /dev/null
+++ b/csv_creation_for_images.py
@@ -0,0 +1,609 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Sat Jul  8 21:15:33 2023
+
+@author: bradenlimb
+"""
+
+#%% Import Modules
+from IPython import get_ipython
+get_ipython().magic('reset -sf')
+import pandas as pd
+from sklearn.preprocessing import normalize
+import numpy as np
+import sys
+import os
+
+import plotly.express as px
+from urllib.request import urlopen
+import json
+from plotly.offline import plot
+
+import geopandas as gpd
+
+import shapefile as shp  # Requires the pyshp package
+import matplotlib.pyplot as plt
+import matplotlib.ticker as ticker
+
+import pickle
+from tqdm import tqdm
+
+import datetime
+begin_time = datetime.datetime.now()
+# asdfa
+
+#%%
+
+# folder_use = '_2023-07-31c Final Results - No Pacific Islands'
+# folder_use = '2023-08-07 Corrected Buffer Life P Level 2'
+# carbon_tax = 20 # $/tonne
+# conc_N = 8 # mg/L
+# conc_P = 1 # mg/L
+
+
+input_data = {}
+input_data['Level 2'] = {
+    'conc_N': 8, # mg/L
+    'conc_P': 1, # mg/L
+    'gray_tech_use': 'Level 2-1, A2O facility'
+    }
+input_data['Level 3'] = {
+    'conc_N': 6, # mg/L
+    'conc_P': 0.2, # mg/L
+    'gray_tech_use': 'Level 4-2, MBR facility'
+    }
+input_data['Level 4'] = {
+    'conc_N': 3, # mg/L
+    'conc_P': 0.1, # mg/L
+    'gray_tech_use': 'Level 4-2, MBR facility'
+    }
+input_data['Level 5'] = {
+    'conc_N': 2, # mg/L
+    'conc_P': 0.02, # mg/L
+    'gray_tech_use': 'Level 5-2, MBR/RO facility'
+    }
+
+## TODO Change Level
+level_use = 5
+folder_use = f'2024-01-03 Level {level_use}'
+conc_N = input_data[f'Level {level_use}']['conc_N'] # mg/L
+conc_P = input_data[f'Level {level_use}']['conc_P'] # mg/L
+gray_tech_use = input_data[f'Level {level_use}']['gray_tech_use']
+
+carbon_tax = 20 # $/tonne
+
+#%% Figure 1
+huc_use = 'HUC6' # HUC6, HUC8, HUC10, or HUC12
+# nutrient = 'Nitrogen'
+# nutrient = 'Phosphorus'
+# with open(f'results/{folder_use}/results_all_point_source_8concN_1concP_postGeoChanges_NorPconc_Land{huc_use}-Percent_GreenCombos_GrayLimits_ConcLimits.pkl', 'rb') as handle:
+with open(f'results/{folder_use}/results_all_level-{level_use}_{huc_use}.pkl', 'rb') as handle:
+    dict_in = pickle.load(handle)
+df_data_in = dict_in['df_data_keep']
+df_data_in.rename(columns = {'HUC 12 Code': 'HUC12'}, inplace=True)
+cols_keep = [huc_use, 'Total Load (kg/yr)', 'Total Annual Flow (L/yr)', 'Mean Conc (mg/L)']
+huc_list = list(df_data_in[huc_use].unique())
+huc_list.remove(huc_list[0])
+huc_list.sort()
+df_data = pd.DataFrame(index = huc_list, columns = cols_keep)
+df_data[huc_use] = df_data.index
+
+# Define how to aggregate various fields
+agg_functions = { 
+                 'Total Load (kg/yr)': 'sum', 
+                 'Total Annual Flow (L/yr)': 'sum',          
+                 }
+
+for nutrient in ['Nitrogen','Phosphorus']:
+    # Create new DataFrame by combining rows with same id values
+    df_data_in_nutrient = df_data_in.loc[df_data_in['Nutrient Type'] == nutrient]
+    df_data = df_data_in_nutrient.groupby([huc_use]).aggregate(agg_functions)
+    df_data['Mean Conc (mg/L)'] = df_data['Total Load (kg/yr)']*1e6/df_data['Total Annual Flow (L/yr)']
+    # df_data[huc_use] = df_data.index
+    
+    # # Save the DataFrame to a CSV file
+    df_data.to_csv(f'results/{folder_use}/nutrient_data_{huc_use}_{nutrient}.csv')
+
+# asdf
+
+#%% Figure 2
+
+land_HUC = huc_use
+
+# Import important HUC12 Data
+with open('inputs/EPA EnviroAtlas/HUC12_Important_Data_v3_withAreasFixed.pkl', 'rb') as handle:
+   df_importantData = pickle.load(handle)
+
+# Remove characters from HUC12 to get other huc values
+df_importantData['HUC6'] = df_importantData['HUC12'].str[:-6]
+df_importantData['HUC8'] = df_importantData['HUC12'].str[:-4]
+df_importantData['HUC10'] = df_importantData['HUC12'].str[:-2]
+
+# df_importantData.loc[df_importantData['area_ha']==0,'area_ha'] = 36*258.999 # Add mean value to each huc12 if area is zero. Source: https://www.epa.gov/rps/step-2-decide-which-watershed-units-screen
+
+# Define a lambda function to compute the weighted mean:
+wm = lambda x: np.average(x, weights=df_importantData.loc[x.index, "area_ha"])
+
+# Define how to aggregate various fields
+agg_functions = { 
+                 'area_ha': 'sum', 
+                 'Cropland_pct': wm,       
+                 'Pasture_pct': wm, 
+                 'Wetlands_pct': wm, 
+                 'StreamLength_km': 'sum', 
+                 'StreamLength_TotalImpaired_km': 'sum', 
+                 'StreamLength_NutImpaired_km': 'sum', 
+                 'StreamLength_TempImpaired_km': 'sum', 
+                 'Wastewater_AnnualFlow_mgy': 'sum', 
+                 'TileDrain_Sum_m2': 'sum', 
+                 'AgBuffer_pct': wm, 
+                 'AgWetlandsPotential_pct': wm, 
+                 'Fertilizer_kgN_ha_yr': wm,
+                 'Elec_$perkWh': wm,
+                 'NoTill_pct': wm,
+                 'CoverCrop_pct': wm,
+                 'N_kg/ha': wm,
+                 'P_kg/ha': wm,
+                 }
+
+# Save original DF as raw DF
+df_importantDataRaw = df_importantData.copy(deep=True)
+
+# Create new DataFrame by combining rows with same id values
+df_importantData = df_importantData.groupby([land_HUC]).aggregate(agg_functions)
+df_importantData.reset_index(inplace=True)
+
+df_importantData.loc[df_importantData['TileDrain_Sum_m2']>0, 'tile_drains'] = 1
+df_importantData.loc[df_importantData['AgBuffer_pct']>0, 'buffer'] = 1
+df_importantData[['Wetlands_pct','AgWetlandsPotential_pct']] = df_importantData[['Wetlands_pct','AgWetlandsPotential_pct']].fillna(0)
+df_importantData.loc[(df_importantData['Wetlands_pct']>0) | (df_importantData['AgWetlandsPotential_pct']>0), 'wetlands'] = 1
+df_importantData.loc[df_importantData['Fertilizer_kgN_ha_yr']>0, 'fertilizer'] = 1
+
+df_importantData[['tile_drains', 'buffer', 'wetlands', 'fertilizer']] = df_importantData[['tile_drains', 'buffer', 'wetlands', 'fertilizer']].fillna(0)
+
+# Save the DataFrame to a CSV file
+df_importantData.to_csv(f'results/{folder_use}/important_data_{land_HUC}.csv')
+
+# sadf
+
+#%% Figure 3
+# Mod 4, slide 88
+# Mod 5, slie 95
+# Mod 5, slie 104
+
+# Import Green Treatment Inputs
+# filename_inputs = 'inputs/model_inputs.xlsx'
+# sheet_green = 'Green Water Treatment Illinois'
+# # sheet_green = 'Green Water Treatment IlliTest'
+# df_green = pd.read_excel(filename_inputs,
+#                                 sheet_name = sheet_green,
+#                                 index_col=1,
+#                                 )
+
+# # Import Results
+# # xslx_path = f'results/{folder_use}/results_summaries_point_source_8concN_1concP_postGeoChanges_NorPconc_Land{land_HUC}-Percent_GreenCombos_GrayLimits_ConcLimits.xlsx'
+# xslx_path = f'results/{folder_use}/results_summaries_level-{level_use}_{huc_use}.xlsx'
+# df_results = pd.read_excel(xslx_path,
+#                         # sheet_name='ghgAGLUC_per_gal-gCO2eq_per_gal',
+#                         # index_col='Land Type'
+#                         )
+
+# Rename the Optimum Names
+
+## Decided to just do this in Excel manually
+
+#%% Figure 5
+
+# Import Green Treatment Inputs
+filename_inputs = 'inputs/model_inputs.xlsx'
+sheet_green = 'Green Water Treatment Illinois'
+# sheet_green = 'Green Water Treatment IlliTest'
+sheet_green = 'Green Water Treatment Illin0'
+df_green = pd.read_excel(filename_inputs,
+                                sheet_name = sheet_green,
+                                index_col=1,
+                                )
+
+green_options = list(df_green['Treatment Name'])[1:]
+
+land_HUC = huc_use
+
+# file_in = f'results/{folder_use}/results_all_point_source_8concN_1concP_postGeoChanges_NorPconc_Land{land_HUC}-Percent_GreenCombos_GrayLimits_ConcLimits.pkl'
+file_in = f'results/{folder_use}/results_all_level-{level_use}_{huc_use}.pkl'
+
+all_options = []
+
+# Import important HUC12 Data
+with open(file_in, 'rb') as handle:
+   dict_resultsAll = pickle.load(handle)
+   
+for data_choice in ['minimal_cost','minimal_emissions']:
+    df_resultsUse =    dict_resultsAll['dict_best_green_option'][data_choice]
+    df_resultsUse2 =    dict_resultsAll['dict_green_treatments']['Cover Crop facility']
+    
+    # # Define a lambda function to compute the weighted mean:
+    # wm = lambda x: np.average(x, weights=df_importantData.loc[x.index, "area_ha"])
+    
+    # Define how to aggregate various fields
+    agg_functions = { 
+        'eGrid': 'first',
+        'N Total Load (kg/yr)': 'sum',
+        'N Annual Flow (L/yr)': 'sum',
+        'N Load Treated (kg/yr)': 'sum',
+        'N Treated Cost ($/yr)': 'sum',
+        'N Treated GWP (tonnes-CO2eq/yr)': 'sum',
+        'P Total Load (kg/yr)': 'sum',
+        'P Annual Flow (L/yr)': 'sum',
+        'P Load Treated (kg/yr)': 'sum',
+        'P Treated Cost ($/yr)': 'sum',
+        'P Treated GWP (tonnes-CO2eq/yr)': 'sum',
+        'Max Treated Cost ($/yr)': 'sum',
+        'Max Treated GWP (tonnes-CO2eq/yr)': 'sum',
+                     }
+    
+    # Save original DF as raw DF
+    df_resultsUseRaw = df_resultsUse.copy(deep=True)
+    
+    # Create new DataFrame by combining rows with same id values
+    df_resultsUse = df_resultsUse.groupby([land_HUC]).aggregate(agg_functions)
+    df_resultsUse.reset_index(inplace=True)
+    
+    # def function ():
+        
+    for option_temp in green_options:
+        option_temp_ab = df_green.loc[df_green['Treatment Name'] == option_temp, 'Abrev'].item()
+        df_resultsUse[option_temp_ab] = 'None'
+        
+        
+    multiHUC_N = []
+    multiHUC_P = []
+    # asdfasdf
+    for hucUse in  tqdm(list(df_resultsUse[land_HUC])):
+        df_tempN = df_resultsUseRaw.loc[(df_resultsUseRaw[land_HUC] == hucUse) & (df_resultsUseRaw['N Load Treated (kg/yr)'] > 0)]
+        n_options = list(set(list(df_tempN['N Option'])))
+        if len(n_options) > 1:
+            multiHUC_N.append(hucUse)
+            if n_options[0] in green_options:
+                n_options[0] = df_green.loc[df_green['Treatment Name'] == n_options[0], 'Abrev'].item()
+            df_resultsUse.loc[df_resultsUse[land_HUC] == hucUse, 'N Option'] = n_options[0]
+            for option_N in n_options[0].split('_'):
+                if df_green.loc[df_green['Abrev'] == option_N, 'N Removal'].item() > 0:
+                    df_resultsUse.loc[df_resultsUse[land_HUC] == hucUse, option_N] = 'Nitrogen'
+        elif len(n_options) == 1:
+            if n_options[0] in green_options:
+                n_options[0] = df_green.loc[df_green['Treatment Name'] == n_options[0], 'Abrev'].item()
+            df_resultsUse.loc[df_resultsUse[land_HUC] == hucUse, 'N Option'] = n_options[0]
+            for option_N in n_options[0].split('_'):
+                if df_green.loc[df_green['Abrev'] == option_N, 'N Removal'].item() > 0:
+                    df_resultsUse.loc[df_resultsUse[land_HUC] == hucUse, option_N] = 'Nitrogen'
+        else:
+            df_resultsUse.loc[df_resultsUse[land_HUC] == hucUse, 'N Option'] = 'None'
+        
+        
+        
+        df_tempP = df_resultsUseRaw.loc[(df_resultsUseRaw[land_HUC] == hucUse) & (df_resultsUseRaw['P Load Treated (kg/yr)'] > 0)]
+        p_options = list(set(list(df_tempP['P Option'])))
+        if len(p_options) > 1:
+            multiHUC_P.append(hucUse)
+            if p_options[0] in green_options:
+                p_options[0] = df_green.loc[df_green['Treatment Name'] == p_options[0], 'Abrev'].item()
+            df_resultsUse.loc[df_resultsUse[land_HUC] == hucUse, 'P Option'] = p_options[0]
+            for option_P in p_options[0].split('_'):
+                if df_green.loc[df_green['Abrev'] == option_P, 'P Removal'].item() > 0:
+                    if df_resultsUse.loc[df_resultsUse[land_HUC] == hucUse, option_P].item() == 'Nitrogen':
+                        df_resultsUse.loc[df_resultsUse[land_HUC] == hucUse, option_P] = 'Both'
+                    else:
+                        df_resultsUse.loc[df_resultsUse[land_HUC] == hucUse, option_P] = 'Phosphorus'
+        elif len(p_options) == 1:
+            if p_options[0] in green_options:
+                p_options[0] = df_green.loc[df_green['Treatment Name'] == p_options[0], 'Abrev'].item()
+            df_resultsUse.loc[df_resultsUse[land_HUC] == hucUse, 'P Option'] = p_options[0]
+            for option_P in p_options[0].split('_'):
+                if df_green.loc[df_green['Abrev'] == option_P, 'P Removal'].item() > 0:
+                    if df_resultsUse.loc[df_resultsUse[land_HUC] == hucUse, option_P].item() == 'Nitrogen':
+                        df_resultsUse.loc[df_resultsUse[land_HUC] == hucUse, option_P] = 'Both'
+                    else:
+                        df_resultsUse.loc[df_resultsUse[land_HUC] == hucUse, option_P] = 'Phosphorus'
+        else:
+            df_resultsUse.loc[df_resultsUse[land_HUC] == hucUse, 'P Option'] = 'None'
+            
+
+    # Replace underscore with plus in the columns
+    df_resultsUse['N Option'] = df_resultsUse['N Option'].str.replace('_', '+')
+    df_resultsUse['P Option'] = df_resultsUse['P Option'].str.replace('_', '+')
+    
+    all_options += list(df_resultsUse['N Option'].unique()) + list(df_resultsUse['P Option'].unique())
+    
+    
+    
+    # Save the DataFrame to a CSV file
+    df_resultsUse.to_csv(f'results/{folder_use}/nutrient_choice_{data_choice}_{land_HUC}.csv')
+    
+    
+all_options = list(set(all_options))
+all_options.remove('None')
+all_options.sort()
+
+# df_alloptions = pd.DataFrame(index = all_options)
+
+# for data_choice in ['minimal_cost','minimal_emissions']:
+#     df_resultsUse =    dict_resultsAll['dict_best_green_option'][data_choice]
+#     df_resultsUse2 =    dict_resultsAll['dict_green_treatments']['Cover Crop facility']
+    
+#     # # Define a lambda function to compute the weighted mean:
+#     # wm = lambda x: np.average(x, weights=df_importantData.loc[x.index, "area_ha"])
+    
+#     # Define how to aggregate various fields
+#     agg_functions = { 
+#         'eGrid': 'first',
+#         'N Total Load (kg/yr)': 'sum',
+#         'N Annual Flow (L/yr)': 'sum',
+#         'N Load Treated (kg/yr)': 'sum',
+#         'N Treated Cost ($/yr)': 'sum',
+#         'N Treated GWP (tonnes-CO2eq/yr)': 'sum',
+#         'P Total Load (kg/yr)': 'sum',
+#         'P Annual Flow (L/yr)': 'sum',
+#         'P Load Treated (kg/yr)': 'sum',
+#         'P Treated Cost ($/yr)': 'sum',
+#         'P Treated GWP (tonnes-CO2eq/yr)': 'sum',
+#         'Max Treated Cost ($/yr)': 'sum',
+#         'Max Treated GWP (tonnes-CO2eq/yr)': 'sum',
+#                      }
+    
+#     # Save original DF as raw DF
+#     df_resultsUseRaw = df_resultsUse.copy(deep=True)
+    
+#     # Create new DataFrame by combining rows with same id values
+#     df_resultsUse = df_resultsUse.groupby([land_HUC]).aggregate(agg_functions)
+#     df_resultsUse.reset_index(inplace=True)
+    
+#     for option in all_options:
+#         df_tempNCount = df_resultsUse.loc(df_resultsUse['N Option'])
+#         df_tempPCount = 
+#         df_alloptions.loc(option, f'{data_choice}_N')
+
+# asdf
+
+#%% eGrid Plot
+huc_use_egrid = 'HUC12'
+
+if huc_use_egrid == 'HUC8':
+    egrid_filename = ''
+elif huc_use_egrid == 'HUC12':
+    egrid_filename = 'inputs/eGrid/HUC12_w_eGrid_v3_smalldf.pkl'
+
+land_HUC = huc_use
+
+with open(egrid_filename, 'rb') as handle:
+    egrid_HUC = pickle.load(handle)    
+    
+egrid_HUC['HUC10'] = egrid_HUC['HUC12'].str[:-2]  
+egrid_HUC['HUC6'] = egrid_HUC['HUC12'].str[:-6]  
+# Define how to aggregate various fields
+agg_functions = { 
+                 'eGrid': 'first', 
+                 'States': 'first',       
+                 }
+
+# Create new DataFrame by combining rows with same id values
+egrid_HUC = egrid_HUC.groupby([land_HUC]).aggregate(agg_functions)
+egrid_HUC.to_csv(f'results/{folder_use}/eGrid_HUC_{land_HUC}.csv')  
+
+
+#%% Find Carbon Financing Percentage
+
+land_HUC = huc_use
+
+# carbon_tax = 20 # $/tonne
+
+# file_in = f'results/{folder_use}/results_all_point_source_8concN_1concP_postGeoChanges_NorPconc_Land{land_HUC}-Percent_GreenCombos_GrayLimits_ConcLimits.pkl'
+file_in = f'results/{folder_use}/results_all_level-{level_use}_{huc_use}.pkl'
+
+# Import important HUC12 Data
+with open(file_in, 'rb') as handle:
+   dict_resultsAll = pickle.load(handle)
+   
+df_gray_useA2O = dict_resultsAll['dict_gray_treatments'][gray_tech_use] 
+   
+dict_comp = {}
+
+
+for data_choice in ['minimal_cost','minimal_emissions']:
+    df_green_use =  dict_resultsAll['dict_best_green_option'][data_choice]
+    
+    df_temp = pd.DataFrame()
+    df_temp[[land_HUC,'Gray Cost ($/yr)', 'Gray Emissions (tonnes-CO2eq/yr)']] = df_gray_useA2O[[land_HUC,'Max Treated Cost ($/yr)','Max Treated GWP (tonnes-CO2eq/yr)']]
+    df_temp[['Green Cost ($/yr)', 'Green Emissions (tonnes-CO2eq/yr)']] = df_green_use[['Max Treated Cost ($/yr)','Max Treated GWP (tonnes-CO2eq/yr)']]
+    df_temp['Emissions Saved (tonnes-CO2eq/yr)'] = df_temp['Gray Emissions (tonnes-CO2eq/yr)'] - df_temp['Green Emissions (tonnes-CO2eq/yr)']
+    
+    df_temp[['N Total Load (kg/yr)', 'N Annual Flow (L/yr)', 'N Load Treated (kg/yr)', 'N Current Mean Conc (mg/L)', 'P Total Load (kg/yr)', 'P Annual Flow (L/yr)', 'P Load Treated (kg/yr)', 'P Current Mean Conc (mg/L)']] = df_green_use[['N Total Load (kg/yr)', 'N Annual Flow (L/yr)', 'N Load Treated (kg/yr)', 'N Current Mean Conc (mg/L)', 'P Total Load (kg/yr)', 'P Annual Flow (L/yr)', 'P Load Treated (kg/yr)', 'P Current Mean Conc (mg/L)']]
+    
+    # Custom lambda function to sum 'N Annual Flow (L/yr)' if 'N Current Mean Conc (mg/L)' > 8
+    custom_sum_N = lambda x: x[df_temp['N Current Mean Conc (mg/L)'] > conc_N].sum()
+    custom_sum_P = lambda x: x[df_temp['P Current Mean Conc (mg/L)'] > conc_P].sum()
+    
+    
+    agg_functions = { 
+        'Gray Cost ($/yr)': 'sum',
+        'Gray Emissions (tonnes-CO2eq/yr)': 'sum',
+        'Green Cost ($/yr)': 'sum',
+        'Green Emissions (tonnes-CO2eq/yr)': 'sum',
+        'Emissions Saved (tonnes-CO2eq/yr)': 'sum',
+        'N Total Load (kg/yr)': custom_sum_N,
+        'N Annual Flow (L/yr)': custom_sum_N,
+        'N Load Treated (kg/yr)': custom_sum_N,
+        'P Total Load (kg/yr)': custom_sum_P,
+        'P Annual Flow (L/yr)': custom_sum_P,
+        'P Load Treated (kg/yr)': custom_sum_P,
+                     }
+    
+    # Save original DF as raw DF
+    df_tempRaw = df_temp.copy(deep=True)
+    
+    # Create new DataFrame by combining rows with same id values
+    df_temp = df_temp.groupby([land_HUC]).aggregate(agg_functions)
+    df_temp.reset_index(inplace=True)
+    
+    # Calculated Needed Treated Load
+    df_temp['N Load Treated Needed (kg/yr)'] = df_temp['N Total Load (kg/yr)'] - conc_N * df_temp['N Annual Flow (L/yr)'] / 1e6
+    df_temp['N Treated Percent (%)'] = df_temp['N Load Treated (kg/yr)'] / df_temp['N Load Treated Needed (kg/yr)'] * 100
+    
+    df_temp['P Load Treated Needed (kg/yr)'] = df_temp['P Total Load (kg/yr)'] - conc_P * df_temp['P Annual Flow (L/yr)'] / 1e6
+    df_temp['P Treated Percent (%)'] = df_temp['P Load Treated (kg/yr)'] / df_temp['P Load Treated Needed (kg/yr)'] * 100
+    
+    
+    # Savings from carbon financing
+    df_temp['Carbon Financing ($/yr)'] = df_temp['Emissions Saved (tonnes-CO2eq/yr)'] * carbon_tax
+    df_temp['Percent Costs (%)'] = df_temp['Carbon Financing ($/yr)'] / df_temp['Green Cost ($/yr)'] * 100
+    
+    # Save the DataFrame to a CSV file
+    df_temp.to_csv(f'results/{folder_use}/carbon_financing_{data_choice}_{land_HUC}_{carbon_tax}.csv')
+    dict_comp[data_choice] = df_temp
+    
+    
+#%% Find the number of facilities for treatment and each nutrient
+df_green_use =  dict_resultsAll['dict_best_green_option']['minimal_cost']
+df_treated_both = df_green_use.loc[(df_green_use['N Current Mean Conc (mg/L)'] >= conc_N) & (df_green_use['P Current Mean Conc (mg/L)'] >= conc_P)]
+df_treated_N = df_green_use.loc[(df_green_use['N Current Mean Conc (mg/L)'] >= conc_N) ]
+df_treated_P = df_green_use.loc[(df_green_use['P Current Mean Conc (mg/L)'] >= conc_P)]   
+df_treated_neither = df_green_use.loc[(df_green_use['N Current Mean Conc (mg/L)'] < conc_N) & (df_green_use['P Current Mean Conc (mg/L)'] < conc_P)]
+ 
+
+#%% Calculate Farmer Incentive Payments
+farmer_incentive = 31 * 2.47105 # $/acre converted to $/ha 
+
+different_costs = ['Buffers','Wetlands', 'Bioreactors']
+df_green_use =  dict_resultsAll['dict_best_green_option']['minimal_cost']
+df_green_use_FI = df_green_use.copy(deep=True)
+df_green_use_FI['max_area'] = df_green_use_FI[['N Treatment Area (ha)','P Treatment Area (ha)']].max(axis=1)
+
+for facility in tqdm(df_green_use_FI.index.tolist()):
+    if (df_green_use_FI.loc[facility, 'N Option'] in different_costs) and (df_green_use_FI.loc[facility, 'P Option'] in different_costs):
+        if df_green_use_FI.loc[facility, 'N Option'] == df_green_use_FI.loc[facility, 'P Option']:
+            df_green_use_FI.loc[facility, 'FI_cost'] = df_green_use_FI.loc[facility, 'max_area'] * farmer_incentive * df_green.loc[df_green['Treatment Name'] == df_green_use_FI.loc[facility, 'N Option'], 'Farmer Incentive'].item()
+        else:
+            print(facility, df_green_use_FI.loc[facility, 'N Option'], df_green_use_FI.loc[facility, 'P Option'])
+    else:
+        df_green_use_FI.loc[facility, 'FI_cost'] = df_green_use_FI.loc[facility, 'max_area'] * farmer_incentive
+
+df_green_use_FI['FI_cost_pct'] = df_green_use_FI['FI_cost'] / df_green_use_FI['Max Treated Cost ($/yr)']
+print(df_green_use_FI['FI_cost'].sum())
+
+#%% Carbon Financing V2
+
+land_HUC = huc_use
+
+# file_in = f'results/{folder_use}/results_facility_point_source_8concN_1concP_postGeoChanges_NorPconc_Land{land_HUC}-Percent_GreenCombos_GrayLimits_ConcLimits.xlsx'
+# file_in = f'results/{folder_use}/results_facility_point_source_{int(conc_N)}concN_{int(conc_P)}concP_postGeoChanges_NorPconc_Land{land_HUC}-Percent_GreenCombos_GrayLimits_ConcLimits.xlsx'
+
+# df_greencost = pd.read_excel(file_in,
+#                                 sheet_name = '72',
+#                                 index_col=0,
+#                                 )
+
+# df_grayA2O = pd.read_excel(file_in,
+#                                 sheet_name = '64',
+#                                 index_col=0,
+#                                 )
+
+# file_in = f'results/{folder_use}/results_all_point_source_8concN_1concP_postGeoChanges_NorPconc_Land{land_HUC}-Percent_GreenCombos_GrayLimits_ConcLimits.pkl'
+file_in = f'results/{folder_use}/results_all_level-{level_use}_{huc_use}.pkl'
+
+
+# Import important HUC12 Data
+with open(file_in, 'rb') as handle:
+   dict_resultsAll = pickle.load(handle)
+   
+data_choice2 = 'minimal_cost'
+df_greencost =  dict_resultsAll['dict_best_green_option'][data_choice2]  
+df_grayA2O = dict_resultsAll['dict_gray_treatments'][gray_tech_use] 
+
+
+
+options_n = list(set(list(df_greencost['N Option'])))
+options_n.sort()
+options_p = list(set(list(df_greencost['P Option'])))
+options_p.sort()
+
+options_all = list(set(options_n + options_p))
+
+options_all_dict = {}
+max_green_cost_name = 'Combined'
+options_all_abr = [max_green_cost_name]
+for option in options_all:
+    if option in green_options:
+        options_all_dict[option] = df_green.loc[df_green['Treatment Name'] == option, 'Abrev'].item()
+        options_all_abr.append(df_green.loc[df_green['Treatment Name'] == option, 'Abrev'].item())
+    else:
+        options_all_dict[option] = option
+        options_all_abr.append(option)
+        
+df_greencost['cost_sum'] =  df_greencost['N Treated Cost ($/yr)'] + df_greencost['P Treated Cost ($/yr)']
+df_greencost['n_cost_portion'] =  df_greencost['N Treated Cost ($/yr)']  / df_greencost['cost_sum'] * df_greencost['Max Treated Cost ($/yr)']
+df_greencost['p_cost_portion'] =  df_greencost['P Treated Cost ($/yr)']  / df_greencost['cost_sum'] * df_greencost['Max Treated Cost ($/yr)']
+  
+df_greencost['emissions_sum'] =  df_greencost['N Treated GWP (tonnes-CO2eq/yr)'] + df_greencost['P Treated GWP (tonnes-CO2eq/yr)']
+df_greencost['n_emissions_portion'] =  df_greencost['N Treated GWP (tonnes-CO2eq/yr)']  / df_greencost['emissions_sum'] * df_greencost['Max Treated GWP (tonnes-CO2eq/yr)']
+df_greencost['p_emissions_portion'] =  df_greencost['P Treated GWP (tonnes-CO2eq/yr)']  / df_greencost['emissions_sum'] * df_greencost['Max Treated GWP (tonnes-CO2eq/yr)']
+
+df_grayA2O['cost_sum'] =  df_grayA2O['N Treated Cost ($/yr)'] + df_grayA2O['P Treated Cost ($/yr)']
+df_grayA2O['n_cost_portion'] =  df_grayA2O['N Treated Cost ($/yr)']  / df_grayA2O['cost_sum'] * df_grayA2O['Max Treated Cost ($/yr)']
+df_grayA2O['p_cost_portion'] =  df_grayA2O['P Treated Cost ($/yr)']  / df_grayA2O['cost_sum'] * df_grayA2O['Max Treated Cost ($/yr)']
+  
+df_grayA2O['emissions_sum'] =  df_grayA2O['N Treated GWP (tonnes-CO2eq/yr)'] + df_grayA2O['P Treated GWP (tonnes-CO2eq/yr)']
+df_grayA2O['n_emissions_portion'] =  df_grayA2O['N Treated GWP (tonnes-CO2eq/yr)']  / df_grayA2O['emissions_sum'] * df_grayA2O['Max Treated GWP (tonnes-CO2eq/yr)']
+df_grayA2O['p_emissions_portion'] =  df_grayA2O['P Treated GWP (tonnes-CO2eq/yr)']  / df_grayA2O['emissions_sum'] * df_grayA2O['Max Treated GWP (tonnes-CO2eq/yr)']
+
+  
+        
+df_financing_Results = pd.DataFrame(index = options_all_abr, 
+                                    columns = ['n_count','n_treated','n_cost_green','n_cost_gray','n_gwp_green','n_gwp_gray','n_cost_cf', 'n_cost_w_cf','p_count','p_treated','p_cost_green','p_cost_gray','p_gwp_green','p_gwp_gray','p_cost_cf', 'p_cost_w_cf',])
+
+
+df_financing_Results.loc[max_green_cost_name,'n_count'] = len(df_greencost)
+df_financing_Results.loc[max_green_cost_name,'n_treated'] = df_greencost['N Load Treated (kg/yr)'].sum()
+df_financing_Results.loc[max_green_cost_name,'n_cost_green'] = df_greencost['n_cost_portion'].sum()
+df_financing_Results.loc[max_green_cost_name,'n_cost_gray'] = df_grayA2O['n_cost_portion'].sum()
+df_financing_Results.loc[max_green_cost_name,'n_gwp_green'] = df_greencost['n_emissions_portion'].sum()
+df_financing_Results.loc[max_green_cost_name,'n_gwp_gray'] = df_grayA2O['n_emissions_portion'].sum()
+
+df_financing_Results.loc[max_green_cost_name,'p_count'] = len(df_greencost)
+df_financing_Results.loc[max_green_cost_name,'p_treated'] = df_greencost['P Load Treated (kg/yr)'].sum()
+df_financing_Results.loc[max_green_cost_name,'p_cost_green'] = df_greencost['p_cost_portion'].sum()
+df_financing_Results.loc[max_green_cost_name,'p_cost_gray'] = df_grayA2O['p_cost_portion'].sum()
+df_financing_Results.loc[max_green_cost_name,'p_gwp_green'] = df_greencost['p_emissions_portion'].sum()
+df_financing_Results.loc[max_green_cost_name,'p_gwp_gray'] = df_grayA2O['p_emissions_portion'].sum()
+
+
+for option in options_all:
+    index_n = df_greencost.loc[df_greencost['N Option'] == option].index.tolist()
+    index_p = df_greencost.loc[df_greencost['P Option'] == option].index.tolist()
+    
+    option_out = options_all_dict[option]
+    
+    df_financing_Results.loc[option_out,'n_count'] = len(index_n)
+    df_financing_Results.loc[option_out,'n_treated'] = df_greencost.loc[index_n,'N Load Treated (kg/yr)'].sum()
+    df_financing_Results.loc[option_out,'n_cost_green'] = df_greencost.loc[index_n,'n_cost_portion'].sum()
+    df_financing_Results.loc[option_out,'n_cost_gray'] = df_grayA2O.loc[index_n,'n_cost_portion'].sum()
+    df_financing_Results.loc[option_out,'n_gwp_green'] = df_greencost.loc[index_n,'n_emissions_portion'].sum()
+    df_financing_Results.loc[option_out,'n_gwp_gray'] = df_grayA2O.loc[index_n,'n_emissions_portion'].sum()
+    
+    df_financing_Results.loc[option_out,'p_count'] = len(index_p)
+    df_financing_Results.loc[option_out,'p_treated'] = df_greencost.loc[index_n,'P Load Treated (kg/yr)'].sum()
+    df_financing_Results.loc[option_out,'p_cost_green'] = df_greencost.loc[index_p,'p_cost_portion'].sum()
+    df_financing_Results.loc[option_out,'p_cost_gray'] = df_grayA2O.loc[index_p,'p_cost_portion'].sum()
+    df_financing_Results.loc[option_out,'p_gwp_green'] = df_greencost.loc[index_p,'p_emissions_portion'].sum()
+    df_financing_Results.loc[option_out,'p_gwp_gray'] = df_grayA2O.loc[index_p,'p_emissions_portion'].sum()
+    
+df_financing_Results['n_cost_cf'] = (df_financing_Results['n_gwp_gray'] - df_financing_Results['n_gwp_green']) * carbon_tax
+df_financing_Results['n_cost_w_cf'] = df_financing_Results['n_cost_green'] - df_financing_Results['n_cost_cf']
+
+df_financing_Results['p_cost_cf'] = (df_financing_Results['p_gwp_gray'] - df_financing_Results['p_gwp_green']) * carbon_tax
+df_financing_Results['p_cost_w_cf'] = df_financing_Results['p_cost_green'] - df_financing_Results['p_cost_cf']
+
+df_financing_Results['n_treated_pct'] = df_financing_Results['n_treated'] / df_financing_Results.loc[max_green_cost_name,'n_treated'] * 100
+df_financing_Results['p_treated_pct'] = df_financing_Results['p_treated'] / df_financing_Results.loc[max_green_cost_name,'p_treated'] * 100
+
+
+df_financing_Results.to_csv(f'results/{folder_use}/carbon_financing_v2_{data_choice2}_{land_HUC}.csv')
diff --git a/green_wastewater_treatment_point-source_data-reduction.py b/green_wastewater_treatment_point-source_data-reduction.py
index 1e95d4e..158bf8c 100644
--- a/green_wastewater_treatment_point-source_data-reduction.py
+++ b/green_wastewater_treatment_point-source_data-reduction.py
@@ -53,11 +53,8 @@
 
 
 farmer_incentive = 31 * 2.47105 # $/acre converted to $/ha - https://theoutcomesfund.com/news-release-2021-environmental-outcomes#:~:text=ANKENEY%2C%20IA%20and%20WASHINGTON%2C%20DC,than%20the%20totals%20from%20the
-# farmer_incentive = 10 * 2.47105 # $/acre converted to $/ha 
-# farmer_incentive = 0
 
 wetland_treatment_area_pct = (2+2*2)/100 # Percent required for the wetland according to 2016 Christenson 10 ways - very conservative. 0.5*2% is range with 2X being the maximum buffer size
-# wetland_treatment_area_pct = 9/100
 
 # Choose the HUC to use for the important data
 HUC_use = 'HUC6' #Choose the HUC level to group by: HUC6, HUC8, HUC10, HUC12
@@ -621,6 +618,7 @@ def treat_nutrients(df_data_in_N, df_data_in_P, df_treatment, treat_conc_N, trea
             df_temp['N Treated Cost ($/yr)'] = df_temp['N Treated Cost ($/yr)'] + df_temp['treated_land_cost']
             df_data_out['N Treated Cost ($/yr)'] = df_temp['N Treated Cost ($/yr)']
     
+    elif (df_treatment['Green/Gray'] == 'Green') & (df_treatment['Farmer Incentive'] == 'No'):
         df_data_out['N Treated Cost ($/yr)'] = df_data_out['N Load Treated (kg/yr)'] * df_treatment['2022Cost kgN']
         # Add land rental costs if wetland is used
         if 'Wetland' in df_treatment['Treatment Name']:
@@ -1211,7 +1209,6 @@ def best_green_summaries(df_in, name):
     plt.show()  # Display the plot
 
 #%% Write the Files to Excel Files
-
 print('Saving Excel Files...')
 
 ## TODO: Change Results Folder Name
@@ -1219,13 +1216,21 @@ def best_green_summaries(df_in, name):
 folder_use = f'{todays_date} Level {min_treat_level}'
 conc_str = f'point_source_{int(treat_conc_N)}concN_{int(treat_conc_P)}concP_postGeoChanges_{NP_output_type}_{land_string}_{green_combos_string}_{limit_gray_string}_{concentration_limit_string}'
 
+conc_str = f'level-{min_treat_level}_{HUC_use}'
+
+
 # Check if the results folder exists
 if not os.path.exists(f'results/{folder_use}'):
     os.makedirs(f'results/{folder_use}') # Create the folder
     
 with pd.ExcelWriter(f'results/{folder_use}/results_summaries_{conc_str}.xlsx') as writer_summaries:
-    df_summary_compare.to_excel(writer_summaries, sheet_name='summary')
+    df_summary_compare.to_excel(writer_summaries, 
+                                sheet_name='summary', 
+                                index = False)
         
+    
+
+    
 raw_results_filename = f'results/{folder_use}/results_facility_{conc_str}.xlsx'
 
 with pd.ExcelWriter(raw_results_filename) as writer_raw:
@@ -1234,22 +1239,26 @@ def best_green_summaries(df_in, name):
     for treatment_green in treatments_green:
         count += 1
         dict_keys[str(count)] = f'{treatment_green}'
+        dict_green_treatments[f'{treatment_green} facility'] = dict_green_treatments[f'{treatment_green} facility'].rename_axis('Facility')
         dict_green_treatments[f'{treatment_green} facility'].to_excel(writer_raw, sheet_name=str(count))
         
     for treatment_gray in treatments_gray:
         count += 1
         dict_keys[str(count)] = f'{treatment_gray}'
+        dict_gray_treatments[f'{treatment_gray} facility'] = dict_gray_treatments[f'{treatment_gray} facility'].rename_axis('Facility')
         dict_gray_treatments[f'{treatment_gray} facility'].to_excel(writer_raw, sheet_name=str(count)) 
     
     count += 1
     dict_keys[str(count)] = 'Green Best Minimal Cost'
+    dict_best_green_option['minimal_cost'] = dict_best_green_option['minimal_cost'].rename_axis('Facility')
     dict_best_green_option['minimal_cost'].to_excel(writer_raw, sheet_name=str(count))
     
     count += 1
     dict_keys[str(count)] = 'Green Best Minimal Emissions'
+    dict_best_green_option['minimal_emissions'] = dict_best_green_option['minimal_emissions'].rename_axis('Facility')
     dict_best_green_option['minimal_emissions'].to_excel(writer_raw, sheet_name=str(count))
     
-    df_keys = pd.DataFrame.from_dict(dict_keys, orient='index', columns=['Values'])
+    df_keys = pd.DataFrame.from_dict(dict_keys, orient='index', columns=['Technology']).rename_axis('Sheet Name')
     df_keys.to_excel(writer_raw, sheet_name='keys')
     
 ## Move the keys sheet to the first sheet
diff --git a/inputs/.DS_Store b/inputs/.DS_Store
index 76eebcd..19c0761 100644
Binary files a/inputs/.DS_Store and b/inputs/.DS_Store differ
diff --git a/inputs/model_inputs.xlsx b/inputs/model_inputs.xlsx
index c3d6e36..9cd349b 100644
Binary files a/inputs/model_inputs.xlsx and b/inputs/model_inputs.xlsx differ
diff --git a/results/.DS_Store b/results/.DS_Store
new file mode 100644
index 0000000..4626046
Binary files /dev/null and b/results/.DS_Store differ
diff --git a/results/2024-01-03 Level 2/.DS_Store b/results/2024-01-03 Level 2/.DS_Store
new file mode 100644
index 0000000..5008ddf
Binary files /dev/null and b/results/2024-01-03 Level 2/.DS_Store differ
diff --git a/results/2024-01-03 Level 2/results_summaries_level-2_HUC6.xlsx b/results/2024-01-03 Level 2/results_summaries_level-2_HUC6.xlsx
new file mode 100644
index 0000000..2271627
Binary files /dev/null and b/results/2024-01-03 Level 2/results_summaries_level-2_HUC6.xlsx differ
diff --git a/results/2024-01-03 Level 3/.DS_Store b/results/2024-01-03 Level 3/.DS_Store
new file mode 100644
index 0000000..5008ddf
Binary files /dev/null and b/results/2024-01-03 Level 3/.DS_Store differ
diff --git a/results/2024-01-03 Level 3/results_summaries_level-3_HUC6.xlsx b/results/2024-01-03 Level 3/results_summaries_level-3_HUC6.xlsx
new file mode 100644
index 0000000..caccf8b
Binary files /dev/null and b/results/2024-01-03 Level 3/results_summaries_level-3_HUC6.xlsx differ
diff --git a/results/2024-01-03 Level 4/.DS_Store b/results/2024-01-03 Level 4/.DS_Store
new file mode 100644
index 0000000..5008ddf
Binary files /dev/null and b/results/2024-01-03 Level 4/.DS_Store differ
diff --git a/results/2024-01-03 Level 4/results_summaries_level-4_HUC6.xlsx b/results/2024-01-03 Level 4/results_summaries_level-4_HUC6.xlsx
new file mode 100644
index 0000000..cccb744
Binary files /dev/null and b/results/2024-01-03 Level 4/results_summaries_level-4_HUC6.xlsx differ
diff --git a/results/2024-01-03 Level 5/.DS_Store b/results/2024-01-03 Level 5/.DS_Store
new file mode 100644
index 0000000..5008ddf
Binary files /dev/null and b/results/2024-01-03 Level 5/.DS_Store differ
diff --git a/results/2024-01-03 Level 5/results_summaries_level-5_HUC6.xlsx b/results/2024-01-03 Level 5/results_summaries_level-5_HUC6.xlsx
new file mode 100644
index 0000000..186dfad
Binary files /dev/null and b/results/2024-01-03 Level 5/results_summaries_level-5_HUC6.xlsx differ
diff --git a/results/Data Dictionary.pdf b/results/Data Dictionary.pdf
new file mode 100644
index 0000000..c5595c9
Binary files /dev/null and b/results/Data Dictionary.pdf differ
diff --git a/results/README.md b/results/README.md
new file mode 100644
index 0000000..14094d1
--- /dev/null
+++ b/results/README.md
@@ -0,0 +1,112 @@
+**Data Dictionary for Results**
+
+_Files:_
+- results\_facility\_level-X\_HUC6.xlsx
+- results\_ summaries\_level-X\_HUC6.xlsx
+- summary\_statistics.xlsx
+
+
+CONTEXT:
+The data included in this data dictionary includes the results that were generated for the publication titled _Greening America's Rivers: The Potential of Climate Financing for Nature-Based Water Quality Improvement._ Similar results files will be generated when the python script is run.
+
+The 'results\_facility\_level-X\_HUC6.xlsx' file contains raw results data from the nutrient trading analysis on the facility level for every green and gray treatment method for the specified treatment level X (levels 2, 3, 4, and 5 evaluated as specified by the [EPA](https://www.epa.gov/system/files/documents/2023-06/life-cycle-nutrient-removal.pdf)).
+
+The 'results\_ summaries\_level-X\_HUC6.xlsx' file contains results data from the nutrient trading analysis aggregated on the national level for every green and gray treatment method for the specified treatment level X (levels 2, 3, 4, and 5 evaluated as specified by the [EPA](https://www.epa.gov/system/files/documents/2023-06/life-cycle-nutrient-removal.pdf)).
+
+The 'summary\_statistics.xlsx' file contains summarized results as presented in the paper titled Greening America's Rivers: The Potential of Climate Financing for Nature-Based Water Quality Improvement.
+
+
+
+VARIABLE DESCRIPTIONS:
+
+File: 'results\_facility\_level-X\_HUC6.xlsx'
+
+Sheet: 'keys'
+- **Sheet Name** – Name of the sheet which includes the technology listed in the 'Technology' column
+- **Technology** – Name of the technology used for nutrient treatment
+
+Sheets: 1+
+- **Facility** – National Pollutant Discharge Elimination System (NPDES) permit number for each wastewater treatment facility
+- **HUC6** – Waterbasin Hydrologic Unit Code as defined by the USGS.
+- **HUC8** – Sub-waterbasin Hydrologic Unit Code as defined by the USGS.
+- **HUC10** – Watershed Hydrologic Unit Code as defined by the USGS.
+- **HUC12** – Sub-watershed Hydrologic Unit Code as defined by the USGS.
+- **eGrid** – Code designating which EPA eGRID electricity region that the facility belongs to.
+- **N Total Load (kg/yr)** – Total nitrogen load released by the wastewater treatment facility each year measured in kilograms per year.
+- **N Annual Flow (L/yr)** – Total annual wastewater flow through each facility measured in liters of water per year.
+- **N Current Mean Conc (mg/L)** – Existing mean concentration of nitrogen in the wastewater measured in milligrams of nitrogen per liter.
+- **N Load Treated (kg/yr)** - Total nitrogen load treated by the nutrient treatment technology each year measured in kilograms per year.
+- **N Treatment Area (ha)** – Number of hectares required for nitrogen treatment using green treatment technologies.
+- **N New Mean Conc (mg/L)** - Mean concentration of nitrogen in the wastewater measured in milligrams of nitrogen per liter after treatment.
+- **N Conc Diff (mg/L)** – Difference in mean nitrogen concentration before and after treatment. Measured in milligrams of nitrogen per liter.
+- **N Treated Cost ($/yr)** – Cost of nitrogen treatment annually. Measured in dollars per year.
+- **N Treated GWP (tonnes-CO2eq/yr)** – Global warming potential of nitrogen treatment annually. Measured in tonnes of CO2 equivalent emissions per year.
+- **P Total Load (kg/yr)** – Total phosphorus load released by the wastewater treatment facility each year measured in kilograms per year.
+- **P Annual Flow (L/yr)** – Total annual wastewater flow through each facility measured in liters of water per year.
+- **P Current Mean Conc (mg/L)** - Mean concentration of phosphorus in the wastewater measured in milligrams of phosphorus per liter.
+- **P Load Treated (kg/yr)** - Total phosphorus load treated by the nutrient treatment technology each year measured in kilograms per year.
+- **P Treatment Area (ha)** – Number of hectares required for phosphorus treatment using green treatment technologies.
+- **P New Mean Conc (mg/L)** - Mean concentration of phosphorus in the wastewater measured in milligrams of phosphorus per liter after treatment.
+- **P Conc Diff (mg/L)** – Difference in mean phosphorus concentration before and after treatment. Measured in milligrams of phosphorus per liter.
+- **P Treated Cost ($/yr)** – Cost of phosphorus treatment annually. Measured in dollars per year.
+- **P Treated GWP (tonnes-CO2eq/yr)** – Global warming potential of phosphorus treatment annually. Measured in tonnes of CO2 equivalent emissions per year.
+- **Max Treated Cost ($/yr)** – Maximum cost of either nitrogen or phosphorus treatment annually, whichever is larger. Measured in dollars per year.
+- **Max Treated GWP (tonnes-CO2eq/yr)** – Maximum global warming potential of either nitrogen or phosphorus treatment annually, whichever is larger. Measured in tonnes of CO2 equivalent emissions per year.
+
+File: 'results\_ summaries\_level-X\_HUC6.xlsx'
+
+- **Treatment** – Name of the technology used for nutrient treatment
+- **Treatment Type –** Identifies if the technology was a Green or Gray treatment option
+- **N Total Load (kg/yr)** – Total nitrogen load released each year measured in kilograms per year.
+- **N Annual Flow (L/yr)** – Total annual wastewater flow measured in liters of water per year.
+- **N Current Mean Conc (mg/L)** – Existing mean concentration of nitrogen in the wastewater measured in milligrams of nitrogen per liter.
+- **N Load Treated (kg/yr)** - Total nitrogen load treated by the nutrient treatment technology each year measured in kilograms per year.
+- **N Treatment Area (ha)** – Number of hectares required for nitrogen treatment using green treatment technologies.
+- **N New Mean Conc (mg/L)** - Mean concentration of nitrogen in the wastewater measured in milligrams of nitrogen per liter after treatment.
+- **N Conc Diff (mg/L)** – Difference in mean nitrogen concentration before and after treatment. Measured in milligrams of nitrogen per liter.
+- **N Treated Cost ($/yr)** – Cost of nitrogen treatment annually. Measured in dollars per year.
+- **N Treated GWP (tonnes-CO2eq/yr)** – Global warming potential of nitrogen treatment annually. Measured in tonnes of CO2 equivalent emissions per year.
+- **P Total Load (kg/yr)** – Total phosphorus load released each year measured in kilograms per year.
+- **P Annual Flow (L/yr)** – Total annual wastewater flow measured in liters of water per year.
+- **P Current Mean Conc (mg/L)** - Mean concentration of phosphorus in the wastewater measured in milligrams of phosphorus per liter.
+- **P Load Treated (kg/yr)** - Total phosphorus load treated by the nutrient treatment technology each year measured in kilograms per year.
+- **P Treatment Area (ha)** – Number of hectares required for phosphorus treatment using green treatment technologies.
+- **P New Mean Conc (mg/L)** - Mean concentration of phosphorus in the wastewater measured in milligrams of phosphorus per liter after treatment.
+- **P Conc Diff (mg/L)** – Difference in mean phosphorus concentration before and after treatment. Measured in milligrams of phosphorus per liter.
+- **P Treated Cost ($/yr)** – Cost of phosphorus treatment annually. Measured in dollars per year.
+- **P Treated GWP (tonnes-CO2eq/yr)** – Global warming potential of phosphorus treatment annually. Measured in tonnes of CO2 equivalent emissions per year.
+- **Max Treated Cost ($/yr)** – Maximum cost of either nitrogen or phosphorus treatment annually, whichever is larger. Measured in dollars per year.
+- **Max Treated GWP (tonnes-CO2eq/yr)** – Maximum global warming potential of either nitrogen or phosphorus treatment annually, whichever is larger. Measured in tonnes of CO2 equivalent emissions per year.
+
+File: 'summary\_statistics.xlsx'
+
+Sheet: 'Statistics'
+- **Cells A1-E15** – Results presented in Table 1 in the paper.
+- **Cells G1-K15** – Statistics on the economic feasibility of using green treatment methods instead of gray treatment methods as a percentage of HUC6 waterbasin.
+
+Sheets: Level X (Where 'X' represents a value between 2-5)
+- **HUC6** – Waterbasin Hydrologic Unit Code as defined by the USGS.
+- **Gray Cost ($/yr)** – Total cost of nutrient treatment in the waterbasin using gray treatment methods.
+- **Gray Emissions (tonnes-CO2eq/yr)** – Total emissions of nutrient treatment in the waterbasin using gray treatment methods.
+- **Green Cost ($/yr)** – Total cost of nutrient treatment in the waterbasin using green treatment methods.
+- **Green Emissions (tonnes-CO2eq/yr)** – Total emissions of nutrient treatment in the waterbasin using green treatment methods.
+- **Emissions Saved (tonnes-CO2eq/yr)** – Total emissions saved by using green treatment methods instead of gray treatment methods for nutrient treatment in the waterbasin.
+- **N Total Load (kg/yr)** – Total nitrogen load released by the wastewater treatment facilities each year in each waterbasin measured in kilograms per year.
+- **N Annual Flow (L/yr)** – Total annual wastewater flow through the wastewater treatment facilities in each waterbasin measured in liters of water per year for plants that treated nitrogen.
+- **N Load Treated (kg/yr)** - Total nitrogen load treated by the nutrient treatment technologies each year in each waterbasin measured in kilograms per year.
+- **P Total Load (kg/yr)** – Total phosphorus load released by the wastewater treatment facilities each year in each waterbasin measured in kilograms per year.
+- **P Annual Flow (L/yr)** – Total annual wastewater flow through the wastewater treatment facilities in each waterbasin measured in liters of water per year for plants that treated phosphorus.
+- **P Load Treated (kg/yr)** - Total phosphorus load treated by the nutrient treatment technologies each year in each waterbasin measured in kilograms per year.
+- **N Load Treatment Needed (kg/yr)** – Total nitrogen load that needed to be treated in each waterbasin to meet desired concentration limits measured in kilograms per year.
+- **N Treated Percent (kg/yr)** – Percent of N Load TreatmentNeeded which was treated in each waterbasin.
+- **P Load Treatment Needed (kg/yr)** – Total phosphorus load that needed to be treated in each waterbasin to meet desired concentration limits measured in kilograms per year.
+- **P Treated Percent (kg/yr)** – Percent of P Load TreatmentNeeded which was treated in each waterbasin.
+- **Carbon Financing ($/yr)** – The total carbon financing in each waterbasin assuming a $20/tonne-CO2e carbon credit cost.
+- **Percent Costs (%)** – The percent of the green treatment costs which could be paid for through carbon financing.
+- **Column T** – The difference between the gray and green treatment costs in each waterbasin in dollars per year.
+- **Column U** – The difference between the gray and green treatment costs including carbon financing in each waterbasin in dollars per year.
+- **Cells W1-Y3** – Calculating the total electricity use for nutrient treatment for gray treatment methods.
+
+Sheets: Level X Totals (Where 'X' represents a value between 2-5)
+- **Cost ($/yr)** – Total costs for each treatment method across all waterbasins in the contiguous United States.
+- **Emissions (tonnes-CO2eq/yr)** - Total emissions for each treatment method across all waterbasins in the contiguous United States.
\ No newline at end of file
diff --git a/results/summary_statistics.xlsx b/results/summary_statistics.xlsx
new file mode 100644
index 0000000..267ff66
Binary files /dev/null and b/results/summary_statistics.xlsx differ