# Python script linking SWAT, MODFLOW and MT3D models # for the SOILPROM project Use case 8: Puck Bay # version under development import csv import numpy as np import flopy import matplotlib.pyplot as plt import flopy.utils as fu import flopy.utils.binaryfile as bf import flopy.utils.formattedfile as ff import subprocess # names of the executable files if each model swat_exe_name = 'SWAT2012.exe' modflow_exe_name = 'MODFLOW-NWT_64.exe' mt3dms_exe_name = 'mt3dms5b.exe' subprocess.run([swat_exe_name]) # spatial coordinates of the model domain for MODFLOW (rectangle) xNW = 451730 yNW = 770900 xSE = 473780 ySE = 747950 # model grid has 459 rows and 441 columns # square gridblocks 50x50 m # reading MODFLOW model information ml = flopy.modflow.Modflow.load('WP.nam',verbose=True) my_nlay = ml.dis.nlay my_nrow = ml.dis.nrow my_ncol = ml.dis.ncol my_nper = ml.dis.nper # reading initial condition for MODFLOW # (hydraulic head) from file WP.hds # after the simulation the file will be overwritten # with values calculated for the current day # and can be used as the initial condition for the next day hfile = bf.HeadFile('WP.hds') # get values from the last time step end_time = hfile.get_times()[-1] init_head = hfile.get_data(totim=end_time) # copy existing array IBOUND my_ibound = ml.bas6.ibound.array.copy() # create new array STRT my_strt = init_head # modify BAS package - new initial condition bas6 = flopy.modflow.mfbas.ModflowBas(model=ml,ibound=my_ibound,strt=my_strt) ml.bas6.write_file() # read recharge date from file output.hru generated by SWAT # for a specific day, e.g. 1 Feb 2001 my_day = 1 my_month = 2 my_year = 2001 # number of HRUs in SWAT model (constant) n_hru = 564 # initialize array with values for HRUs swat_data = np.zeros((n_hru,6)) with open ('output.hru') as out_file: out_reader = csv.reader(out_file, delimiter=' ', skipinitialspace=True) for i in range(9): # skip the heading next(out_reader,None) while True: row = next(out_reader,None) if row == None: break r_hru = int(row[1])-1 r_month = int(row[5]) r_day = int(row[6]) r_year = int(row[7]) #print(r_day, r_month, r_year) r_area = float(row[8]) # area of HRU [km2] r_perc = float(row[17]) # recharge [mm/d] r_evap = float(row[20]) # ET from groundwater [mm/d] r_no3 = float(row[64]) # load N-NO3 [kg/ha/d] # check the date if (r_day == my_day) and (r_month == my_month) and (r_year == my_year): swat_data[r_hru,0]=r_area swat_data[r_hru,1]=r_perc swat_data[r_hru,2]=r_evap swat_data[r_hru,3]=r_no3 swat_data[r_hru,4]=0 swat_data[r_hru,5]=0 if r_perc > r_evap: # net recharge [m/h] swat_data[r_hru,4]=(r_perc-r_evap)/24000 # N-NO3 concentration in recharge [mg/dm3] swat_data[r_hru,5]=r_no3/swat_data[r_hru,4]/240 # calculate average values for SWAT model area area_total = np.sum(swat_data[:,0]) rch_aver = np.sum(np.multiply(swat_data[:,0],swat_data[:,4]))/area_total no3_aver = np.sum(np.multiply(swat_data[:,0],swat_data[:,3]))/area_total no3_conc_aver = 0 if rch_aver > 0: no3_conc_aver = no3_aver/rch_aver/240 print('SWAT model area: ',str(area_total),' km2') print('Average recharge: ',str(rch_aver), 'm/h') print('Average N-NO3 conc. in recharge: ',str(no3_conc_aver), 'mg/dm3') # read files mapping SWAT HRUs to MODFLOW grid cells # hru_dhru.txt and grid_dhru.txt # reading file hru_dhru.txt # indices of HRUs corresponding to specific DHRUs idhh = [] with open ('hru_dhru.txt') as out_file: out_reader = csv.reader(out_file, delimiter='\t', skipinitialspace=True) # read number of DHRUs (divides HRU) row = next(out_reader,None) n_dhru = int(row[0]) # read number of HRUs row = next(out_reader,None) #my_hru = int(row[0]) # read heading row = next(out_reader,None) for i in range(n_dhru): row = next(out_reader,None) idhh.append(int(row[2])) grid_data = [] with open ('grid_dhru.txt') as out_file: out_reader = csv.reader(out_file, delimiter='\t', skipinitialspace=True) # read the number of records row = next(out_reader,None) n_grid = int(row[0]) # read the number of DHRUs (divided HRUs) row = next(out_reader,None) #n_dhru = int(row[0]) # read the number of rows row = next(out_reader,None) #my_row = int(row[0]) # read the number of columns row = next(out_reader,None) #my_ncol = int(row[0]) # read heading row = next(out_reader,None) for i in range(n_grid): row = next(out_reader,None) grid_data.append(row) #my_grid = int(grid_data[0][0]) #print(my_grid) #print(ml.dis.get_lrc(my_grid)[0][1]) # preparing array with recharge values # in all grid cells my_mask = ml.bas6.ibound.array[2,:,:].copy() my_rch = np.zeros((my_nrow,my_ncol)) # recharge equal to the average value for the whole area my_rch[my_mask>0] = rch_aver my_no3_load = np.zeros((my_nrow,my_ncol)) # nitrate load equal to the average value for the whole area my_no3_load[my_mask>0] = no3_conc_aver*rch_aver for igrid in grid_data: my_grid = int(igrid[0]) # find the row and column index of the given grid cell my_row = ml.dis.get_lrc(my_grid)[0][1] my_col = ml.dis.get_lrc(my_grid)[0][2] #print(my_grid,my_row,my_col) # check if the cell is active if my_mask[my_row,my_col] > 0: # modify recharge and concentration grid_area = float(igrid[1]) overlap = float(igrid[3]) coeff = overlap/grid_area #print(grid_area,overlap,coeff) # find DHRU index, indexing from 0 my_dhru = int(igrid[2])-1 # find HRU index, indexing from 0 my_hru = idhh[my_dhru]-1 #print(my_dhru,my_hru) hru_rch = swat_data[my_hru,4] hru_cno3 = swat_data[my_hru,5] new_rch = my_rch[my_row,my_col]+(hru_rch-rch_aver)*coeff new_no3_load = (my_no3_load[my_row,my_col]\ +(hru_cno3*hru_rch-rch_aver*no3_conc_aver)*coeff) #print(my_rch[my_row,my_col]) #print(hru_rch,new_rch) my_rch[my_row,my_col] = new_rch my_no3_load[my_row,my_col] = new_no3_load #print(my_rch[my_row,my_col]) #input('Press Enter') my_no3_conc = np.zeros((my_nrow,my_ncol)) my_no3_conc[my_rch>0] = np.divide(my_no3_load[my_rch>0],my_rch[my_rch>0]) # adding new recharge values to the model flopy.modflow.mfrch.ModflowRch(ml, nrchop=3, ipakcb=40, \ rech=my_rch, unitnumber=10) ml.rch.write_file() ml.exe_name = modflow_exe_name # running MODFLOW program success, buff = ml.run_model() print(success) print(buff) # reading MT3D model mt = flopy.mt3d.mt.Mt3dms.load('WP.mt_nam',exe_name=mt3dms_exe_name, verbose=True) # specify MODFLOW file which is the basis of the transport simulation mt.ftlfilename = 'WP.ftl' # read N-NO3 concentration (from the previous step) # file WP_NO3.ucn, which will be later overwritten with the results for the current time step my_c0 = np.zeros((my_nlay, my_nrow, my_ncol)) cfile = bf.UcnFile('WP_NO3.ucn') # control listing of time steps #for t in cfile.get_times(): # print(str(t)) # results from the last time step end_time = cfile.get_times()[-1] #print(str(end_time)) for ilay in range(my_nlay): my_c0[ilay,:,:] = cfile.get_data(totim=end_time,mflay=ilay) #cN_2 = cfile.get_data(totim=end_time,mflay=2) #cN_4 = cfile.get_data(totim=end_time,mflay=4) # writing new initial condition for MT3D (concentrations) to file WP.btn mt.btn.sconc = [] u3d = flopy.utils.util_array.Util3d(mt, (my_nlay, my_nrow, my_ncol), np.float32, my_c0, name='sconc1', locat=mt.btn.unit_number[0], array_free_format=False) mt.btn.sconc.append(u3d) mt.btn.write_file() # writing new values of N-NO3 load to file WP.ssm mt.ssm.crch = [] t2d = flopy.utils.util_array.Transient2d(mt, (my_nrow, my_ncol), np.float32, my_no3_conc, name='crch1', locat=mt.ssm.unit_number[0], array_free_format=False) mt.ssm.crch.append(t2d) mt.ssm.write_file() # modify the first row in WP.ssm, # necessary to read the file in subsequent steps ssm_file = open('WP.ssm','rt') # open file ssm_lines = ssm_file.readlines() # read all lines ssm_file.close() # close the file ssm_lines[0]=' F F T F F T F F F F\n' ssm_file = open('WP.ssm','wt') # open file for writing for line in ssm_lines: # write line ssm_file.write(line) ssm_file.close() # close the file # running MT3D computer program success, buff = mt.run_model() print(success) print(buff) # reading fluxes between grid cells ffile = bf.CellBudgetFile('WP.cbc') #ffile.list_records() # control listing of time steps #for t in ffile.get_times(): # print(str(t)) # results from the last time step end_time = ffile.get_times()[0] #print(str(end_time)) # flow in x firection (W-E) flow1 = ffile.get_data(text='FLOW RIGHT FACE', totim=end_time,full3D=True) # flow in y direction (N-S) flow2 = ffile.get_data(text='FLOW FRONT FACE', totim=end_time,full3D=True) # flow in z direction flow3 = ffile.get_data(text='FLOW LOWER FACE', totim=end_time,full3D=True) # read the calculated concentration of N-NO3 # we are interested in layers "2" and "4" (indexing from 0) # kontrolny wydruk kroków czasowych z pliku wynikowego cfile = bf.UcnFile('WP_NO3.ucn') # control listing of time steps #for t in cfile.get_times(): # print(str(t)) # results from the last time step end_time = cfile.get_times()[0] #print(str(end_time)) cN_2 = cfile.get_data(totim=end_time,mflay=2) cN_4 = cfile.get_data(totim=end_time,mflay=4) # control plots of hydraulic head and concentration hdata = bf.HeadFile('WP.hds') # get data from the last time step end_time = hdata.get_times()[-1] my_heads = hdata.get_data(totim=end_time) #fig = plt.figure(figsize=(10, 6)) #ax = fig.add_subplot(1, 2, 1, aspect="equal") #modelmap = flopy.plot.PlotMapView(model=ml, layer=0, ax=ax) #pa1 = modelmap.plot_array(new_head, masked_values=[-999.], cmap='rainbow') #cb = plt.colorbar(pa1, shrink=0.5) plot_heads = np.copy(my_heads[2, :, :]) plot_heads[my_heads[2, :, :] < -900.0] = np.nan plot1 = plt.figure(1) plt.subplot(1, 1, 1, aspect='equal') plt.title('Head distribution (m)') plt.imshow(plot_heads, cmap='rainbow', vmin = 0.0) plt.colorbar() #contours = plt.contour(np.flipud(head[0, :, :]), levels=levels, extent=extent, zorder=10) #plt.clabel(contours, inline=1, fontsize=10, fmt='%d', zorder=11) #plt.show() plot2 = plt.figure(2) cdata = fu.UcnFile('WP_NO3.UCN') end_time = cdata.get_times()[-1] conc = cdata.get_data(totim=end_time, mflay=2) conc[my_heads[2, :, :] < -900.0] = np.nan #conc.plot(totim=end_time, colorbar='NO3 concentration (mg/l)', cmap='Blues', vmin=0.) plt.imshow(conc, cmap='rainbow', vmin = 0.0) plt.title('Concentration distribution (mg/l)') plt.colorbar() plt.show()