Matlab model Christ and Goltz

Matlab model for the calculation of capture zones of wells by Christ and Goltz

% CaptureZone2000.m
% Script for calculating the capture zone of the wells 
% Akacievej, Floeng. 
%
% Based on the paper of Christ and Goltz, J. Hydrology, v262, 224-244, 2002.
%
% Code based on a script written by MSc students
% Lise Celine Pedersen and Kristine Duelund Vilsgaard 
%
% Philip Binning, August 2013.
% Klaus Mosthaf, August 2014.
clear all
close all
clc

%% Variable definitions

K=1e-3; % Conductivity [m/s], lowest one: 8e-5 m/s
GradH=-0.0008; % Hydraulic head gradient
U=-K*GradH; % Darcy flux [m/s]
Ho=20.0; % Reference head at K11 [m]
B=25; % Aquifer thickness [m] 15-35m
alpha=25/360*2*pi; % 2*3.14; % Angle of flow with x-axis in radians, N-NE


Q(1)=23000/365/24/3600;%2.5/3600; % Pump rate at Hedehusene vestre [m3/s]
Q(2)=80000/365/24/3600;%12/3600; % Pump rate at Hedehusene østre [m3/s]
Q(3)=15/3600; % Pump rate Akacievej (remedial pump) [m3/s]
Q(4)=8500/365/24/3600;%/3600; % Pump rate in 106390 Marbjerg VV [m3/s]
Q(5)=70000/365/24/3600;%/3600; % Pump rate in Position of Fløng WW 1 [m3/s]
Q(6)=3500000/365/24/3600;%/3600; % Pump rate in Solhøj WW [m3/s], 3.7-6.5 Mio. m3/yr
Q(7)=0/365/24/3600;%/3600; % Pump rate at Nyfløng [m3/s]
Q(8)=830000/365/24/3600;%115/3600; % Pump rate in Katrinebjerg [m3/s], 0.8-1.2 Mio. m3/yr
Q(9)=10000/365/24/3600;%3600; % Pump rate in Soderup [m3/s], 0.8-1.2 Mio. m3/yr
Q(10)=4200000/365/24/3600; % Pump rate in Marbjerg kildeplads [m3/s]
Q(11)=70000/365/24/3600; % Pump rate in Fløng 2 [m3/s]
Q(12)=13500/365/24/3600; % Pump rate in 106400 [m3/s]

% for i=1:8
%     Q(i)=0; % Set all pumping rates to zero
% end

% Coordinate system and grid
xmin=-3440;
xmax=7710;
xstep=10;
ymin=-3010;
ymax=4370;
ystep=10;
ixo=round(-xmin/xstep);  % index of point at (x,y)=(0,0)
jyo=round(-ymin/ystep);

X=(xmin:xstep:xmax);
Y=(ymin:ystep:ymax);
[X,Y]=meshgrid(X,Y);

% Positions of wells 

x(1)=2060; y(1)=-35; % Position of Hedehusene vestre vv (K11)
x(2)=2525; y(2)=300; % Position of Position of Hedehusene østre vv
x(3)=2520-1000; y(3)=1985-1000; % Position of remedial pump
x(4)=1305-1000; y(4)=2598-1000; % Position of 106390 (Marbjerg)
x(5)=2300-1000; y(5)=2742-1000; % Position of Fløng WW
x(6)=4913-1000; y(6)=-1100-1000; % Position of Solhøj WW
x(7)=1790; y(7)=550; % Position of 106400
x(8)=5682-1000; y(8)=7040-1000; % Position of origin
x(9)=1925; y(9)=3200; % Position of remedial pump
x(10)=-1200; y(10)=1950; % Position of Marbjerg Vandværk
x(11)=1860; y(11)=1960; % Position of remedial pump
x(12)=5910; y(12)=480; % Position of remedial pump
x(13)=2552-1000; y(13)=1956-1000; % Position of remedial pump
namew{1}='Hedehusene vestre';
notew{1}='22.000 m^3/yr';
namew{2}='Hedehusene østre';
notew{2}='100.000 m^3/yr';
namew{3}='Akacievej';
notew{3}='22.000 m^3/yr';
namew{4}='Marbjerg By'; %(106390)
notew{4}='10.000 m^3/yr';
namew{5}='Fløng 1'; %(106387)
notew{5}='140000 m^3/yr';
namew{6}='Solhøj WW';
notew{6}='3.700.000 m^3/yr';
namew{7}='Nyfløng';
notew{7}='0 m^3/yr';
namew{8}='Katrinebjerg';
notew{8}='1.000.000 m^3/yr';
namew{9}='Sodereup';
notew{9}='10.000 m^3/yr';
namew{10}='Marbjerg';
notew{10}='4.800.000 m^3/yr';
namew{11}='Fløng 2'; %(106387)
notew{11}='70000 m^3/yr';
namew{12}='Høje Thorstrup';
notew{12}='105.000 m^3/yr';
namew{13}='Akac_out';
notew{13}='105.000 m^3/yr';


% Definition of the complex field W
Z=X+1i*Y;
W=-U*Z*exp(-1i*alpha);
for j=1:12 % loop over number of wells
    zj=x(j)+1i*y(j);
    W = W + Q(j)/(2*pi*B)*log(Z-zj);
end

% Set constant C in Christ equation (1) so that Head at (0,0) = Ho.
% Remember here that H=phi/K and that real(w)=phi.

W = W + Ho*K - real(W(ixo,jyo)); 

%% Background map
scrsz=get(0,'ScreenSize');
figure('Units', 'normalized', 'OuterPosition',[0.01, 0.01, .8, .8])

%I=imread('FloengHedehusene_noAnnotations.jpg');
%I=imread('DWplants_cut.png');
I=imread('VandvaerkLarge_cut.png');
imagesc([xmin xmax], [ymax ymin], I)
set(gca,'YDir','normal')

%% Plot contours and streamlines
figure(1)
hold on
[C1,h1]=contour(X,Y,imag(W),45,'blue');  
[C2,h2]=contour(X,Y,real(W)/K,'color',[0,0.7,0],'linewidth',2); % Divide by K, since it is velocity potential
set(h2,'LevelStep',0.5);
set(h2,'ShowText','on','TextStep',get(h2,'LevelStep')*2);
axis equal;
axis([xmin xmax ymin ymax]);
legend('show','Streamlines','Potentials','Location','NorthWest');
xlabel('Distance [m]'); ylabel('Distance [m]');

% Plot wells and contaminated site.
for j=1:12
    plot(x(j),y(j),'.r', 'MarkerSize', 20)
    text(x(j)+50,y(j)+90,namew{j}, 'FontSize', 16)
%    text(x(j)+50,y(j)-80,notew{j}, 'FontSize', 14)
end

Back to Useful tools