TITLE '2D solute flow across Lichtaard, NW Friesland, Netherlands '
SELECT
stages=10
contours=10
VARIABLES     { system variables }
pressure
omega
DEFINITIONS   { parameter definitions }
g=vector(0,10)
gy=10
k=1.2e-10
mu=1e-3
D=9.5e-10
rho_0=1000
n=0.35
rho=rho_0*exp(0.69*omega)			{density; function of omega (salt concentration)}
J=-rho*D*grad(omega)			{diffusive flux}
q=-(k/mu)*(grad(pressure)+rho*g)
w=2.14e-3*erf(y/(8/3))			{concentration; function of y}
p={stated in attached word document}
fac=10**((stage-10)/3)
EQUATIONS     { PDE's, one for each variable }
pressure:div(rho*q*fac)=0
omega:div(rho*omega*q*fac)+div(J)=0
BOUNDARIES    { The domain definition }
  REGION 1 "total area"
    START(0,0)
	natural(pressure)=p {Here I need an integral expression}
	natural(omega)=w
	line to (0,-32)
	natural(pressure)=0 {The base boundary is closed; no groundwaterflow through it}
	value(omega)=
	line to (1250,-32)
	natural(pressure)=p {Here I need an integral expression}
	natural(omega)=w
	line to (1250,0)
	value(pressure)=0
	value(omega)=
	line to (1150,0)
	value(pressure)=p
	value(omega)=w
	line to (1150,-1.5)
	value(pressure)=0
	value(omega)=
	line to (900,-1.5)
	value(pressure)=p
	value(omega)=w
	line to (900,0)
	value(pressure)=0
	value(omega)=
	line to (650,0)
	value(pressure)=p
	value(omega)=w
	line to (650,-1.5)
	value(pressure)=0
	value(omega)=
	line to (450,-1.5)
	value(pressure)=p
	value(omega)=w
	line to (450,0)
	value(pressure)=0
	value(omega)=
	line to (300,0)
	value(pressure)=p
	value(omega)=w
	line to (300,-1.5)
	value(pressure)=0
	value(omega)=
	line to (100,-1.5)
	value(pressure)=p
	value(omega)=w
	line to (100,0)
	value(pressure)=0
	value(omega)=
	line to finish
 REGION 2 "deklaag, layer 1"
  k=2.3e-14
  n=0.10
  START(0,0)

	line to (0,-2.0)

	line to (1250,-2.0)

	line to (1250,0)

	line to (1150,0)

	line to (1150,-1.5)

	line to (900,-1.5)

	line to (900,0)

	line to (650,0)

	line to (650,-1.5)

	line to (450,-1.5)

	line to (450,0)

	line to (300,0)

	line to (300,-1.5)

	line to (100,-1.5)

	line to (100,0)

	line to finish
 REGION 3 "fijnzandig wvp,layer 2"
 k=1.2e-13
 n=0.20
  START(0,-2.0)

	line to (0,-9.0)

	line to (1250,-9.0)

	line to (1250,-2.0)

	line to finish
 REGION 4 "grofzandig wvp,layer 3"
 k=1.2e-10
 n=0.35
  START(0,-9.0)

	line to (0,-32)

	line to (1250,-32)

	line to (1250,-9.0)

	line to finish
REGION 5 "bodem afgraving 1"
 k=7e-14
 n=0.35
  START(100,-1.5)

	line to (100,-2.0)

	line to (300,-2.0)

	line to (300,-1.5)

	line to finish

REGION 6 "bodem afgraving 2"
 k=7e-14
 n=0.35
  START(450,-1.5)

	line to (450,-2.0)

	line to (650,-2.0)

	line to (650,-1.5)

	line to finish

REGION 7 "bodem afgraving 3"
 k=7e-14
 n=0.35
  START(900,-1.5)

	line to (900,-2.0)

	line to (1150,-2.0)

	line to (1150,-1.5)

	line to finish


PLOTS


elevation(pressure) FROM (0,0) to (0,-32)
elevation(pressure) from (1250,0) to (1250,-32)
contour(pressure) zoom (300,0,350,-5) painted
contour(omega) zoom (300,0,350,-5) painted


END