#########################################################################################
# Introduce heterogeneity 
g.region rows=50 cols=100 res3=50 t=500 b=0 n=1000 s=0 w=0 e=2000 -p3

r.mapcalc expression="mask = if(col()> 20 && col() < 60 && row()>20 && row()<35, 1, 2)" --overwrite

r.mapcalc expression="phead = if(col() == 1, 30, 29)" --overwrite
r.mapcalc expression="status = if(col() == 1  || col() == 100, 2, 1)" --overwrite
r.mapcalc expression="top_unconf = 100.0" --overwrite
r.mapcalc expression="null = 0.0" --overwrite
r.mapcalc expression="poros = 0.15" --overwrite

r.mapcalc expression="hydcond = if (mask == 1, 0.0001, 0.0025)" --overwrite

r.gwflow top=top_unconf bottom=null phead=phead status=status hc_x=hydcond hc_y=hydcond s=poros output=gwresult type=unconfined dt=864000000 vx=velX vy=velY error=0.000000000000001 maxit=100000000000 --overwrite

r.mapcalc "magnitude = sqrt(velX^2 + velY^2)" --overwrite
r.mapcalc "direction = atan(velX, velY)" --overwrite

d.rast.arrow map=direction color=black magnitude_map=magnitude skip=20 grid=none scale=50

###########################################################################################
# Confined aquifer
g.region rows=50 cols=100 res3=50 t=500 b=0 n=1000 s=0 w=0 e=2000 -p3

r.mapcalc expression="phead = if(col() == 1, 30, 29)" --overwrite
r.mapcalc expression="status = if(col() == 1  || col() == 100, 2, 1)" --overwrite
r.mapcalc expression="hydcond = 0.0025" --overwrite
r.mapcalc expression="storativ = 0.00001" --overwrite
r.mapcalc expression="top_conf = 10" --overwrite
r.mapcalc expression="null = 0.0" --overwrite

r.gwflow top=top_conf bottom=null phead=phead status=status hc_x=hydcond hc_y=hydcond s=storativ output=gwresult type=confined dt=864000000 vx=velX vy=velY error=0.000000000000001 maxit=100000000000 --overwrite