# -*- coding: utf-8 -*- """ Created on Mon Sep 30 12:26:17 2019 @author: jaza0001 Compare the numerical simulation and analytical solution of partially penetrating river freshwater lens """ import os import sys import numpy as np import flopy.utils.binaryfile as bf import matplotlib.pyplot as plt try: import flopy except: fpth=os.path.abspath(os.path.join('..','..')) import flopy print (sys.version) print('numpy version:{}'.format(np.__version__)) print('flopy version:{}'.format(flopy.__version__)) workspace=os.path.join('test A1') if not os.path.exists(workspace): os.makedirs(workspace) print(workspace) # SEAWAT Model Class modelname='Partial' swt=flopy.seawat.Seawat(modelname,exe_name='swt_v4x64',model_ws=workspace) # Input parameters ************************************************************ # Experiment A-1--------------------------------------------------------------- eta_r = 0.06478 # river penetration depth into the sand (m) eta_a = 0.24534 # depth of the sand tank base below the riverbed (m) W_r = 0.05035 # half of the river width (m) B_r = 0.01041 # screen (resistive material or clogging layer) thickness (m) eta_sb = 0.32506 # saltwater thickness at the saltwater reservoir located at x_b (m) eta_sL = eta_r+B_r+eta_a # saltwater thickness the lens tip (m) x_b = 1.04499 # horizontal distance between the inlet and outlet reservoirs (m) rho_f = 995.99970 # freshwater density (kg/m3) rho_s = 1029.99141 # saltwater density (kg/m3) K_fa = 228.47847 # freshwater hydraulic conductivity of the sand (to be used in SEAWAT) (m/d) k=0.41779 # factor accounting the area of the sand cross-sectional area occupied by the screen K_fc = K_fa*k # clogging layer (or screen) freshwater hydraulic conductivity (to be used in SEAWAT) (m/d) K_sa = K_fa*rho_s/rho_f # saltwater hydraulic conductivity of the sand (to be used in analytical solution) (m/d) K_sc = K_fc*rho_s/rho_f # clogging layer (or screen) saltwater hydraulic conductivity (to be used in analytical solution) (m/d) # # Experiment A-2--------------------------------------------------------------- # eta_r = 0.11523 # eta_a = 0.19511 # W_r = 0.05035 # B_r = 0.01041 # eta_sb = 0.32517 # eta_sL = eta_r+B_r+eta_a # x_b = 1.04499 # rho_f = 995.99105 # rho_s = 1030.04512 # K_fa = 234.02688 # k=0.41779 # K_fc = K_fa*k # K_sa = K_fa*rho_s/rho_f # K_sc = K_fc*rho_s/rho_f # # Experiment A-3--------------------------------------------------------------- # eta_r = 0.16523 # eta_a = 0.14510 # W_r = 0.05035 # B_r = 0.01041 # eta_sb = 0.32495 # eta_sL = eta_r+B_r+eta_a # x_b = 1.04499 # rho_f = 996.00004 # rho_s = 1030.03677 # K_fa = 237.74182 # k=0.41779 # K_fc = K_fa*k # K_sa = K_fa*rho_s/rho_f # K_sc = K_fc*rho_s/rho_f # # Experiment B-1--------------------------------------------------------------- # eta_r = 0.06494 # eta_a = 0.24507 # W_r = 0.09934 # B_r = 0.01088 # eta_sb = 0.32452 # eta_sL = eta_r+B_r+eta_a # x_b = 0.99517 # rho_f = 996.00334 # rho_s = 1029.99467 # K_fa = 222.84316 # k=0.39154 # K_fc = K_fa*k # K_sa = K_fa*rho_s/rho_f # K_sc = K_fc*rho_s/rho_f # # Experiment B-2--------------------------------------------------------------- # eta_r = 0.11533 # eta_a = 0.19495 # W_r = 0.09934 # B_r = 0.01088 # eta_sb = 0.32482 # eta_sL = eta_r+B_r+eta_a # x_b = 0.99517 # rho_f = 995.99935 # rho_s = 1029.99973 # K_fa = 222.37869 # k=0.39154 # K_fc = K_fa*k # K_sa = K_fa*rho_s/rho_f # K_sc = K_fc*rho_s/rho_f # # Experiment B-3--------------------------------------------------------------- # eta_r = 0.16531 # eta_a = 0.14504 # W_r = 0.09934 # B_r = 0.01088 # eta_sb = 0.32482 # eta_sL = eta_r+B_r+eta_a # x_b = 0.99517 # rho_f = 995.99944 # rho_s = 1030.00047 # K_fa = 212.64131 # k=0.39154 # K_fc = K_fa*k # K_sa = K_fa*rho_s/rho_f # K_sc = K_fc*rho_s/rho_f # ============================================================================= # SEAWAT simulation # ============================================================================= # MODFLOW Discretization Packag------------------------------------------------ Lx=W_r+B_r+x_b Lz=eta_sb delv=0.002 delr=0.002 delc=1 nlay=np.int(Lz/delv) ncol=np.int(Lx/delr) nrow=1 nper=1 top=Lz botm=np.linspace(top-delv,0,nlay) perlen=2 nstp=48 tsmult=1.02 dis=flopy.modflow.ModflowDis(swt, nlay=nlay, nrow=nrow, ncol=ncol, nper=nper, delr=delr, delc=delc, laycbd=0, top=top, botm=botm, perlen=perlen, nstp=nstp, tsmult=tsmult, steady=False, itmuni=4, lenuni=2, extension='dis', unitnumber=None, filenames=None, xul=None, yul=None, rotation=None, proj4_str=None, start_datetime=None) # MODFLOW Basic Package-------------------------------------------------------- hdry=-999 hf=eta_r+B_r+eta_a hfl=nlay-np.int(hf/delv+0.5) if hfl <0: hfl = 0 rbl=np.int((eta_r+B_r)/delv+0.5)+hfl rbc=np.int((W_r+B_r)/delr+0.5) nlc=np.int(B_r/delv+0.5) ncc=np.int(B_r/delr+0.5) ibound=np.ones((nlay,nrow,ncol),dtype=np.int32) ibound[:rbl,:,:rbc]=0 ibound[:,0,-1]=-1 strt=np.ones((nlay,nrow,ncol),dtype=np.float32)*eta_sb bas=flopy.modflow.ModflowBas(swt, ibound=ibound, strt=strt, ifrefm=True, ixsec=False, ichflg=False, stoper=None, hnoflo=hdry, extension='bas', unitnumber=None, filenames=None) # MODFLOW Layer Property Flow Package------------------------------------------ laytyp=1 hk=np.ones((nlay,nrow,ncol),dtype=np.float32)*K_fa # aquifer hydraulic conductivity-in SEAWAT freshwater hydraulic conductivity should be used hk[rbl-nlc:rbl,:,rbc-ncc:rbc]=K_fc # clogging layer hydraulic conductivity-in SEAWAT freshwater hydraulic conductivity should be used ss=1e-06 sy=0.24 lpf=flopy.modflow.ModflowLpf(swt, laytyp=laytyp, layavg=0, chani=1.0, layvka=0, laywet=0, ipakcb=53, hdry=hdry, iwdflg=0, wetfct=0.1, iwetit=1, ihdwet=0, hk=hk, hani=1.0, vka=hk, ss=ss, sy=sy, vkcb=0.0, wetdry=-0.01, storagecoefficient=False, constantcv=False, thickstrt=False,nocvcorrection=False, novfc=False, extension='lpf',unitnumber=None, filenames=None) # MODFLOW General-Head Boundary Package (for riverbank)------------------------ K_ev=K_fc*K_fa*(B_r+delv*0.5)/(K_fa*B_r+K_fc*delv*0.5) # equivalent K for vertical flow from riverbed K_eh=K_fc*K_fa*(B_r+delr*0.5)/(K_fa*B_r+K_fc*delr*0.5) # equivalent K for horizontal flow from riverbank cond_rbnk=delc*delv*K_eh/(B_r+delr*0.5) # conductance of riverbank cond_rbd=delc*delr*K_ev/(B_r+delv*0.5) # conductance of riverbed ghba_dtype=np.dtype([('k', '0.2 and d_p<=0.5: a1,a2 = [0.538,-0.387] else: a1,a2 = [0.538,-0.387] elif W_pn>1 and W_pn<=3.5: if d_p<=0.2: a1,a2 = [0.819,-1.340] elif d_p>0.2 and d_p<=0.5: a1,a2 = [0.672,-0.542] elif d_p>0.5 and d_p<=0.9: a1,a2 = [0.567,-0.330] kapa = np.exp(-np.pi*W_pn) Gama_flat = 1/(2*(1+(1/np.pi)*np.log(2/(1-np.sqrt(kapa))))) # conductance for flat river (no penetration) Gama_p = Gama_flat*(1+a1*d_p+a2*d_p**2) # conductance with penetarated river Gama_c = Gama_p/(1+2*(B_r/W_p)*(K_sa/K_sc)*Gama_p) # conductance with clogging layer # Assume x_B < x_L , x_L freshwater lens extent-------------------------------- a = -1/(2*K_sa*Gama_c**2*(1-rho_f/rho_s)) b = (x_b-x_B)+(B_r+eta_a)/Gama_c c = 0.5*K_sa*(eta_sb**2-rho_f/rho_s*eta_sL**2-(1-rho_f/rho_s)*(B_r+eta_a)**2) q_s1 = (-b+np.sqrt(b**2-4*a*c))/(2*a) q_s2 = (-b-np.sqrt(b**2-4*a*c))/(2*a) if q_s1<0: # only negative q_s is accepted as negative sign indicates flow towards river q_s=q_s1 elif q_s2<0: q_s=q_s2 x_L = 0.5*K_sa*(eta_sb**2-eta_sL**2)/q_s+x_b eta_sB = B_r+eta_a-q_s/(K_sa*Gama_c*(1-rho_f/rho_s)) # if assume x_B > x_L , x_L freshwater lens extent----------------------------- if x_L0.2 and d_p<=0.5: a1,a2 = [0.538,-0.387] elif W_pn>1 and W_pn<=3: if d_p<=0.2: a1,a2 = [0.819,-1.340] elif d_p>0.2 and d_p<=0.5: a1,a2 = [0.672,-0.542] elif d_p>0.5 and d_p<=0.9: a1,a2 = [0.567,-0.330] kapa = np.exp(-np.pi*W_pn) Gama_flat = 1/(2*(1+(1/np.pi)*np.log(2/(1-np.sqrt(kapa))))) # conductance for flat river (no penetration) Gama_p = Gama_flat*(1+a1*d_p+a2*d_p**2) # conductance with penetarated river Gama_c = Gama_p/(1+2*(B_r/W_p)*(K_sa/K_sc)*Gama_p) # conductance with clogging layer # Assume x_B < x_L , x_L freshwater lens extent-------------------------------- a = -1/(2*K_sa*Gama_c**2*(1-rho_f/rho_s)) b = (x_b-x_B)+(B_r+eta_a)/Gama_c c = 0.5*K_sa*(eta_sb**2-rho_f/rho_s*eta_sL**2-(1-rho_f/rho_s)*(B_r+eta_a)**2) q_s1 = (-b+np.sqrt(b**2-4*a*c))/(2*a) q_s2 = (-b-np.sqrt(b**2-4*a*c))/(2*a) if q_s1<0: q_s=q_s1 elif q_s2<0: q_s=q_s2 x_L = 0.5*K_sa*(eta_sb**2-eta_sL**2)/q_s+x_b eta_sB = B_r+eta_a-q_s/(K_sa*Gama_c*(1-rho_f/rho_s)) # if assume x_B > x_L , x_L freshwater lens extent----------------------------- if x_L