University of California
Groundwater

# Tutorial: 1D ADE transport simulation

## Demo of 1D transport equation

` `

### Problem definition

The domain of the problem will be 5 km long.

```L = 5000;
```

We will discretize the domain into 50 m linear line elements. Therefore the coordinates of the points are

```p = ( 0:50:L )';
Np = size(p, 1);
```

Next we define the 1D mesh. To generate a uniform 1D mesh we do not need any special software, however we need to create a structure variable, to hold the mesh information, similar to the one we use in more complex meshes. The first 2 lines of code are not used and descibe the 0 dimension elements e.g. boundary points. (You could set MSH(1,1) = []; instead)

```MSH(1,1).elem(1,1).type = 'Bndpnt';
MSH(1,1).elem(1,1).id = [1;Np];
MSH(2,1).elem(1,1).type = 'line';
MSH(2,1).elem(1,1).id = [(1:Np-1)' (2:Np)'];
Nel = size(MSH(2,1).elem(1,1).id,1);
```

For initial conditions we will use a constant concentration of of 50 mg/L at the first node with id 1

```CC = [1 50];
```

The initial distribution of the concentration will be zero

```Cinit = zeros(Np, 1);
```

and we will run the simulation for 150 years with yearly step

```T = (0:150)'*365;
```

Since we do not have any water flows we set a vector of zeros

```F = sparse(Np,1);
```

and we will use crank Niclolson scheme

```wmega = 0.5;
```

Last we will define the few simulation options

```opt.dim=1; % This is the dimension of the problem
opt.el_type='line'; %the element type
opt.el_order='linear';% the element order (linear is the only valid option)
opt.assemblemode='vect';% theis is the mode. (use always vectorized option)
opt.capacmode='consistent';% option regarding the capacitance matrix (other option is 'lumped')
```

### Constant transport parameters, no decay, no retardation

First we will illustrate the most simple transport case with all transport properties constant. In addition we will assume no retardation and decay

```aL = 500; %[m] longitudinal dispersivity
v = 0.15; %[m/day] velocity
lambda = 0; %[1/day] radioactive decay
K_d = 0; %[m^3/Kg] equilibrium distribution coefficient
rho_b = 1;% bulk density
Dm = 1.1578e-004;%[m^2/day] Molecular diffusion coefficient
theta = ones(Nel,1);
```

To assemble the mass and dispersivity matrix we call the function

```[Dglo Mglo c]= Assemble_LHS_std(p, MSH(2,1).elem(1,1).id,...
aL, v, rho_b, K_d, lambda, theta, Dm, CC, opt);
```

... and we are ready to solve the transport equation

```C1 = SteadyFlowTransport(Mglo, Dglo, F, Cinit, T, c, wmega);
```

Plotting concentration profiles and the breakthrough curve at the outlet in matlab is easy

```figure('Position',[100 100 660 220])
subplot(1,2,1);
surf(p/1000,T(2:end)/365,C1,'edgecolor','none')
view(0,90)
axis([0 5 0 150])
xlabel('Distance [km]')
ylabel('Time [years]')
colorbar
subplot(1,2,2)
plot(T(2:end)/365,C1(:,end))
xlabel('Time [years]')
ylabel('Concentration at the outlet [mg/L]')

```

### Constant transport parameters, with decay and retardation

Now lets examine the effect of the radioactive decay on the previous problem. lambda is defined as ln(2)/half-life. We will iterate through various lambda coefficients which correspond to different half-life times.

```half_life = [10:10:100]*365;
lambda = log(2)./half_life;
clear lgnd
lgnd{1,1} = 'no decay';
for i = 1:length(lambda)
[Dglo Mglo c]= Assemble_LHS_std(p, MSH(2,1).elem(1,1).id,...
aL, v, rho_b, K_d, lambda(i), theta, Dm, CC, opt);
c_temp = SteadyFlowTransport(Mglo, Dglo, F, Cinit, T, c, wmega);
C_l(:,i) = c_temp(:,end);
lgnd{i+1,1} = [num2str(half_life(i)/365) ' yr half-life'];
end
figure('Position',[100 100 560 420])
plot(T(2:end)/365,C1(:,end),'--r')
hold on
plot(T(2:end)/365,C_l)
xlabel('Time [years]')
ylabel('Concentration at the outlet [mg/L]')
title('Decay effect')
legend(lgnd,'Location','NorthWest')
hold off
```

Next we will set the the decay equal to zero and solve for various retardation values

```lambda = 0;
K_d = [0.5:0.5:4];
clear lgnd
lgnd{1,1} = 'R = 1';
for i = 1:length(K_d)
[Dglo Mglo c]= Assemble_LHS_std(p, MSH(2,1).elem(1,1).id,...
aL, v, rho_b, K_d(i), lambda, theta, Dm, CC, opt);
c_temp = SteadyFlowTransport(Mglo, Dglo, F, Cinit, T, c, wmega);
C_r(:,i) = c_temp(:,end);
lgnd{i+1,1} = ['R = ' num2str(1+K_d(i))];
end
clf
plot(T(2:end)/365,C1(:,end),'--r')
hold on
plot(T(2:end)/365,C_r)
xlabel('Time [years]')
ylabel('Concentration at the outlet [mg/L]')
title('Retardation effect')
legend(lgnd,'Location','NorthWest')
```

### Variable transport parameters

Last we will see how to use variable parameters. Let the velocity be a function of x as follows:

```V_fnc = inline('0.5*exp(-((x-2500)/500).^2)+0.15');
```

The velocity profile along the domain is defined as

```V_nd = V_fnc(p);
figure('Position',[100 100 300 300])
plot(p, V_fnc(p));
```

We solve this by simply passing the velocity matrix as agrument

```K_d = 0;
[Dglo Mglo c]= Assemble_LHS_std(p, MSH(2,1).elem(1,1).id,...
aL, V_fnc(p), rho_b, K_d, lambda, theta, Dm, CC, opt);
C_v = SteadyFlowTransport(Mglo, Dglo, F, Cinit, T, c, wmega);
figure('Position',[100 100 560 420])
subplot(2,2,1);
surf(p/1000,T(2:end)/365,C1,'edgecolor','none')
view(0,90)
axis([0 5 0 150])
xlabel('Distance [km]')
ylabel('Time [years]')
title('Constant parameters')
colorbar
subplot(2,2,2);
surf(p/1000,T(2:end)/365,C_v,'edgecolor','none')
view(0,90)
axis([0 5 0 150])
xlabel('Distance [km]')
ylabel('Time [years]')
title('Variable Velocity')
colorbar
subplot(2,2,3);
plot(T(2:end)/365,C1(:,end))
xlabel('Time [years]')
ylabel('Concentration at the outlet [mg/L]')
axis([0 150 0 50])
subplot(2,2,4);
plot(T(2:end)/365,C_v(:,end))
xlabel('Time [years]')
ylabel('Concetration at the outlet [mg/L]')
axis([0 150 0 50])
```

Now lets define some other transport properties as function of x. Note that typically these do not make sence in real cases.

The longitudinal dispersivity will be a linear function of x which varies from 200 m at the inlet to 5000 m at the outlet.

```aL_fnc = inline('200+x*0.96');
```

We solve this in similar way

```[Dglo Mglo c]= Assemble_LHS_std(p, MSH(2,1).elem(1,1).id,...
aL_fnc(p), V_fnc(p), rho_b, K_d, lambda, theta, Dm, CC, opt);
C_aL = SteadyFlowTransport(Mglo, Dglo, F, Cinit, T, c, wmega);
```

Now lets define a variable decay constant. We assume that the decay is zero for the first half of the domain and then the decay constant corresponds to 30 years half-life.

```lmd_fnc = inline('(x>2500)*log(2)/(30*365)');
[Dglo Mglo c]= Assemble_LHS_std(p, MSH(2,1).elem(1,1).id,...
aL_fnc(p), V_fnc(p), rho_b, K_d, lmd_fnc(p), theta, Dm, CC, opt);
C_lmd = SteadyFlowTransport(Mglo, Dglo, F, Cinit, T, c, wmega);
figure('Position',[100 100 560 420])
subplot(2,2,1);
surf(p/1000,T(2:end)/365,C_aL,'edgecolor','none')
view(0,90)
axis([0 5 0 150])
xlabel('Distance [km]')
ylabel('Time [years]')
title('... and Variable Dispersivity')
colorbar
subplot(2,2,2);
surf(p/1000,T(2:end)/365,C_lmd,'edgecolor','none')
view(0,90)
axis([0 5 0 150])
xlabel('Distance [km]')
ylabel('Time [years]')
title('... and Variable Decay')
colorbar
subplot(2,2,3);
plot(T(2:end)/365,C_aL(:,end))
xlabel('Time [years]')
ylabel('Concentration at the outlet [mg/L]')
axis([0 150 0 50])
subplot(2,2,4);
plot(T(2:end)/365,C_lmd(:,end))
xlabel('Time [years]')
ylabel('Concentration at the outlet [mg/L]')
axis([0 150 0 50])
```

Note that the code to solve the transport is just few lines. The majority of the code above is about plotting

` `
Webmaster Email: thharter@ucdavis.edu