The following guide shows how to prepare the files for a steady and transient state simulation. Since this is only a demo for the transient simulation we will use a very small domain. Note also to make things less complicated we assume that the groundwater aquifer is bounded on the top by the ground surface. This is obviously not true and one should consider a steady state representation of the water table instead for top.

In addition, the following shows the input file preparation when a single processor is used. For multi processor simulations check this page

The C2VSim model can be downloaded from the DWR repository. The download includes only the inputs and one need to run the model to create the outputs.

By default, the velocity is not printed. To print the velocity field you need to set a filename for the VELOUTFL parameter of the main Groundwater input file. (For version 1 and 1.01 the name of this file is C2VSimFG_Groundwater1974.dat). For this tutorial we assume that we set the name of the velocity output as C2VSimFG_GW_VELOUTFL.out.

There are also two parameters that modify how the velocity is printed. The FACTVROU is a multiplier to the velocity that was set by default to 0.000022957. This has to be changed to FACTVROU=1000.0. The second parameter UNITVROU it is used only for displaying the correct units and does not affect the numeric values. Since we change the FACTVROU=1000.0 for consistency reasons set UNITVROU = 10^3FT/MON.

You may notice that the default value was AC-FT/MON e.g. L^3/T but we changed to a unit that is L/T. After personal communication I confirmed that the velocity is printed in L/T not in L^3/T. If we set the FACTVROU = 1 then the velocity is printed in FT/MON. However, there are parts of the flow domain where the velocity is very low and the printing accuracy does not allow to capture these small ranges. For that reason we set the FACTVROU = 1000 which increases the printed accuracy. If this value is increased to a very large multiplier value e.g. 10^6 then a few of the velocity numbers become very large which exceed the allowable printing width.

Before getting into the input files of Ichnos there are two output data from C2VSim that are time-consuming to read. These are the Groundwater heads and the Velocity file. It is recommended to read those in advance and save them as Matlab variables which makes the loading significantly faster.

Set up the paths

c2vsim_path = fullfile('path','to','C2VsimV1','c2vsim-working');
input_files_path = fullfile('path','to','input','files');

Read the data

Heads = readIWFM_headalloutput(fullfile(c2vsim_path,'Results',"C2VSimFG_GW_HeadAll.out"), size(ND,1), 4, 505, 1);
VelOut = readIWFM_Velocity(fullfile(c2vsim_path,'Results','C2VSimFG_VETOUTFL.out'), 32537, 4, 1000);

Save the data

save(fullfile(c2vsim_path,'GW_Heads.mat'), Heads);
save(fullfile(c2vsim_path,'GW_Velocity.mat'), VelOut);

For the domain files we have converted the nodes and elements files to a polygon shapefile. In a GIS software, the elements shapefile was dissolved to a single polygon that is the outline of the domain. It has also been re-projected to 3310.

First we read the shapefile and split the polygon because the outline consists of several rings.

c2vsim_outline = shaperead('C2VsimMesh_Outline_3310.shp'); 
[Xs, Ys] = polysplit(c2vsim_outline.X, c2vsim_outline.Y);

Next we write the outline to the file after each subpolygon has been simplified. The simplify threshold requires some experimentation.

simplify_threshold = 1000;
fid = fopen(fullfile(input_files_path, 'c2vsimFG_outline.ich') ,'w');
for ii = 1:length(Xs)
    [xx,yy] = reducem(Xs{ii,1}', Ys{ii,1}', simplify_threshold);
    fprintf(fid, '%d %d\n',[length(xx) ispolycw(xx,yy)]);
    fprintf(fid, '%.3f %.3f\n', [xx yy]');
end
fclose(fid);

Since the top and bottom elevations are both defined on the same XY nodes we will write one file containing both elevations.

In the following snippets we use functions from gwtools

Read the Nodes where the elevation is set.

ND = readIWFM_Nodes(fullfile(c2vsim_path,'Preprocessor','C2VSimFG_Nodes.dat'));

The coordinates of the original file are in the EPSG:26910. Here we convert them to EPSG:3310.

[lat,lon] = projinv(projcrs(26910),ND(:,1), ND(:,2));
[X_3310, Y_3310] = projfwd(projcrs(3310),lat, lon);

Read the Elevation data. Note that the elevation data are in feet, so we convert them to meters

Strat = readIWFM_Stratigraphy(fullfile(c2vsim_path,'Preprocessor','C2VSimFG_Stratigraphy.dat'), size(ND,1), 4, 105);
TopElev = Strat(:,2) * 0.3048; 
Bottom_elev = (Strat(:,2) - sum(Strat(:,3:end),2)) * 0.3048;

Write both top and bottom elevation data in one file

fid = fopen(fullfile(input_files_path,'c2vsim_TopBot_CLOUD.ich'),'w');
fprintf(fid, 'CLOUD\n');
fprintf(fid, '%.2f %.2f\n', [4000, 3]);
fprintf(fid, '%.3f %.3f %.4f %.4f\n', [X_3310 Y_3310 TopElev Bottom_elev]');
fclose(fid);

Processor polygons are not needed for single processing simulation

To read the velocity output file, we use one of the Matlab functions of gwtools. Note that the last argument of the function is the FACTVROU multiplier. The function will read and divide the velocity to convert it to FT/MON.

VelOut = readIWFM_Velocity(fullfile(c2vsim_path,'Results','C2VSimFG_GW_VELOUTFL.out'), 32537, 4, 1000);

The output is a structure that contains the coordinates of the 2D element barycenters where the velocity is defined. However, the coordinates are in ft as well as in a different projection system.

[lat,lon] = projinv(projcrs(26910),VelOut.ND(:,1)*0.3048, VelOut.ND(:,2)*0.3048);
[BCX_3310, BCY_3310] = projfwd(projcrs(3310),lat, lon);

First we will run a steady state particle tracking. We will average the velocities of the last decade of the simulation (10/2005 - 9/2015).

for ii = 1:size(VelOut.VX,3)
    Vx_av(:,ii) = mean(VelOut.VX(:,385:504,ii),2);
    Vy_av(:,ii) = mean(VelOut.VY(:,385:504,ii),2);
    Vz_av(:,ii) = mean(VelOut.VZ(:,385:504,ii),2);
end

Calculate the elevation of the velocity nodes.

First compute the elevations of the element mesh nodes

ND_ELEV(:,1) = Strat(:,2)*0.3048;
ND_ELEV(:,2) = ND_ELEV(:,1) - sum(Strat(:,3:4),2)*0.3048;
ND_ELEV(:,3) = ND_ELEV(:,2) - sum(Strat(:,5:6),2)*0.3048;
ND_ELEV(:,4) = ND_ELEV(:,3) - sum(Strat(:,7:8),2)*0.3048;
ND_ELEV(:,5) = ND_ELEV(:,4) - sum(Strat(:,9:10),2)*0.3048;

Then calculate the elevation of the element barycenters, the element diameter and the ratio between xy/z element extend and the average element thickness

msh = readIWFM_Elements(fullfile(c2vsim_path,'Preprocessor','C2VSimFG_Elements.dat'), 32537, 142);
for ii = 1:size(msh,1)
    nv = 4;
    if msh(ii,5) == 0
        nv = 3;
    end
    for j = 1:4
        t = mean(ND_ELEV(msh(ii,2:nv+1),j));
        b = mean(ND_ELEV(msh(ii,2:nv+1),j+1));
        MSH_ELEV(ii,j) = (t + b)/2;
        av_thick(ii,j) =  mean(ND_ELEV(msh(ii,2:nv+1),j) - ND_ELEV(msh(ii,2:nv+1),j+1));
    end
    DIAM(ii,1) = max(pdist([X_3310(msh(ii,2:nv+1)) Y_3310(msh(ii,2:nv+1))]));
end

Convert the velocity units from FT/MON to m/day.

Vx_av = Vx_av .* (0.3048/30.437);
Vy_av = Vy_av .* (0.3048/30.437);
Vz_av = Vz_av .* (0.3048/30.437);

The velocity values, especially the z components are very small e.g. 10^-5. Therefore, we will multiply the values with a very large constant multiplier. During particle tracking we have to reverse this by setting the multiplier equal to 1/mult.

Prepare and reshape the data for printing:

nvel_p = size(MSH_ELEV,1)*size(MSH_ELEV,2);
mult = 1000000;
VEL_DATA = [...
    repmat(BCX_3310,4,1) ... % X
    repmat(BCY_3310,4,1) ... % X
    reshape(MSH_ELEV, nvel_p, 1) ... % Z
    zeros(nvel_p,1) ... % Proc ID
    repmat(DIAM,4,1) ... % DIAM 
    reshape(bsxfun(@rdivide, DIAM, av_thick), nvel_p, 1) ... % RATIO
    reshape(Vx_av, nvel_p, 1)*mult ... % VX
    reshape(Vy_av, nvel_p, 1)*mult ... % VY
    reshape(Vz_av, nvel_p, 1)*mult ... % VZ
];

Last we can write the data as

• ASCII

fid = fopen(fullfile(input_files_path,'c2vsim_VelCLOUD_SS_000.ich'),'w');
fprintf(fid, '%.3f %.3f %.3f %d %.2f %.2f %.5f %.5f %.5f\n', VEL_DATA');
fclose(fid);

Using ASCII file has the advantage that we can easily examine the printed files. Yet the loading times are quite large for very large velocity fields.

• HDF5

HDF5 format reduces the loading times and we can print the data in an efficient HDF5 format as follows:

fname = fullfile(input_files_path,'c2vsim_VelCLOUD_SS_000.h5');
h5create(fname,'/XYZDR',[size(VEL_DATA,1) 5], 'Datatype','single');
h5write(fname, '/XYZDR', VEL_DATA(:,[1 2 3 5 6]));
h5create(fname,'/PROC',[size(VEL_DATA,1) 1], 'Datatype','uint32');
h5write(fname, '/PROC', VEL_DATA(:,4));
h5create(fname,'/VXYZ',[size(VEL_DATA,1) 3], 'Datatype','single');
h5write(fname, '/VXYZ', VEL_DATA(:,[7 8 9]));

The following is optional but recommended. It will create a file that shows the connectivity of the cloud. Then every time step the algorithm will search for the nearest point and use only the connected points to the nearest point for interpolation.

Let's create the 2D connection graph of the mesh. For each element we will make a list of the elements that that is connected to

msh_v = msh(:,2:5); % Isolate the mesh vertices
clear msh_graph
msh_graph{size(msh,1),1} = [];
for ii = 1:size(msh,1)
    nd = msh_v(ii,:);
    nd(:, nd == 0) = [];
    con_elem = [];
    for jj = 1:length(nd)
        [II, JJ] = find(msh_v == nd(jj));
        con_elem = [con_elem;II];
    end
    con_elem = unique(con_elem);
    con_elem(con_elem == ii,:) = [];
    msh_graph{ii,1} = con_elem;
end

Each velocity node in the 3D cloud depends on the nodes of the elements of the same layer, the layer above and the layer below.

clear mesh3D_graph
nel = size(msh_graph,1);
mesh3D_graph{nel,4} = [];
for ii = 1:size(msh,1)
    for jj = 1:4
        if jj == 1
            mesh3D_graph{ii,jj} = [jj*nel+ii; (jj-1)*nel+msh_graph{ii}; jj*nel+msh_graph{ii}];
        elseif jj == 4
            mesh3D_graph{ii,jj} = [(jj-2)*nel+ii; (jj-2)*nel+msh_graph{ii}; (jj-1)*nel+msh_graph{ii}];
        else
            mesh3D_graph{ii,jj} = [(jj-2)*nel+ii; jj*nel+ii; ...
                (jj-2)*nel+msh_graph{ii}; (jj-1)*nel+msh_graph{ii}; jj*nel+msh_graph{ii}];
        end
    end
end
mesh3D_graph = reshape(mesh3D_graph, nvel_p, 1);

Print the file. (see here for details)

fid = fopen(fullfile(input_files_path, 'c2vsim_SS_05_15_000.grph'), 'w');
for ii = 1:length(mesh3D_graph)
    fprintf(fid, '%.2f %.2f %.2f', VEL_DATA(ii,1:3));
    fprintf(fid, ' %d %d', [length(mesh3D_graph{ii}) 1]);
    fprintf(fid, ' %d', [mesh3D_graph{ii}' ii]-1);
    fprintf(fid, '\n');
end
fclose(fid);

Next we will generate a few particles to trace within the steady state field. Typically the locations of the particles depend on the study. For the purpose of this tutorial we will manually select 20 locations in x and y

plot(c2vsim_outline.X, c2vsim_outline.Y)
particles_locations = ginput(20);

Calculate the top elevation at these locations

Ftop = scatteredInterpolant(X_3310, Y_3310, TopElev, 'linear');
pTop = Ftop(particles_locations(:,1),particles_locations(:,2));

In each location we will release 4 particles with 10 m apart in the vertical direction. As you can see the particles are grouped in 4 entities with ids 10, 20, 30, 40. Within each entity the streamline ids span from 1-20.

particleData = [ ...
    10*ones(length(pTop),1) (1:length(pTop))' particles_locations pTop - 10];
particleData = [particleData; ...
    20*ones(length(pTop),1) (1:length(pTop))' particles_locations pTop - 20];
particleData = [particleData; ...
    30*ones(length(pTop),1) (1:length(pTop))' particles_locations pTop - 30];
particleData = [particleData; ...
    40*ones(length(pTop),1) (1:length(pTop))' particles_locations pTop - 40];

Print the particles to a file.

fid = fopen(fullfile(input_files_path,'CV_particles.ich'),'w');
fprintf(fid, '# Test particle file\n');
fprintf(fid, '#\n');
fprintf(fid, '%d %d %.3f %.3f %.3f\n', particleData');

In addition to single particles we will generate particles from hypothetical wells.

We will generate 3 wells with Eid {100, 200, 300} and top 10 m below the top function and 50 m screen length. Note that we specify the release time to 8022 days. This will be ignored in the steady state.

plot(c2vsim_outline.X, c2vsim_outline.Y)
well_locations = ginput(3);
wtop = Ftop(well_locations(:,1), well_locations(:,2));
well_data = [ (100:100:300)' well_locations wtop - 10 wtop - 60 8022*ones(length(wtop),1)];

Write the data to file. We will release 20 particles per well organized into 4 layers at a distance of 10 m from the x,y well location

fid = fopen(fullfile(input_files_path,'CV_wells.ich'),'w');
fprintf(fid, '# Test well file\n');
fprintf(fid, '#\n');
fprintf(fid, '# Cloud point tutorial\n');
fprintf(fid, '# Npart Nlay Rad\n');
fprintf(fid, '%d %d %.2f\n', [25 5 10.0]); 
fprintf(fid, '# Eid X Y T B RT\n');
fprintf(fid, '%d %.3f %.3f %.3f %.3f %.1f\n', well_data');
fclose(fid);

To run the code set up the two configuration files, the main configuration file c2vsim_CLOUD_SS_config.ini and the velocity configuration file c2vsim_CLOUD_SS_vf.ini.

For the step configuration we deactivate the StepSizeTime and the nStepsTime by setting them to 0 so that we can control the step size without taking into account the field velocity speed.

⚠️
Make sure that the output folder exists before running

Then run the code as:

ichnos.exe -c c2vsim_CLOUD_SS_config.ini

Alternatively you can run the simulation tracking using more processors. Make sure that the you set the options RunAsThread = 1 and nThreads = 1

C:\"Program Files\Microsoft MPI"\Bin\mpiexec.exe -n 4 ichnos.exe -c c2vsim_CLOUD_SS_config.ini

To read the results in matlab you can use one of the functions of gwtools

S = readICHNOStraj('CLOUD_SS_iter_0000.traj');

The following snippet plots the streamlines.

clf
hold on
for ii = 1:length(S)
    plot3(S(ii,1).p(:,1), S(ii,1).p(:,2), S(ii,1).p(:,3),'-')
end
plot(c2vsim_outline.X, c2vsim_outline.Y, 'linewidth',2)
axis equal
axis off

We run the code with 2 step configurations. One with StepSize = 50 and one with StepSize = 500

If we zoom into one of the wells we can see the cyan streamlines that correspond to 50 m step and the yellow streamlines(500 m). It is obvious that the particle trajectories with large Euler steps appear rather wiggly.

In the next figure we focus on one of the locations where we release 4 streamlines at different depths. The cyan and yellow streamlines correspond to Euler with very small and large step respectively. The purple and the red correponds to RK2 and RK4 using the same large step, while the blue is based on RK45. We can see that the RK45 results in similar streamlines as the very fine Euler

In the next sections we will demonstrate how to prepare the input files for transient velocity fields. For the purpose of this demo we will use the flow field of the last decade of simulation.

Se above how to read the velocity field. The units of the velocity of the VelOut variable are FT/month and the following snippet converts them to m/day

days_per_month = eomday(VelOut.YMD(:,1), VelOut.YMD(:,2));
Vx_m_day = VelOut.VX;
Vy_m_day = VelOut.VY;
Vz_m_day = VelOut.VZ;
for ii = 1:length(days_per_month)
    Vx_m_day(:,ii,:) = 0.3048 .* Vx_m_day(:,ii,:) ./ days_per_month(ii);
    Vy_m_day(:,ii,:) = 0.3048 .* Vy_m_day(:,ii,:) ./ days_per_month(ii);
    Vz_m_day(:,ii,:) = 0.3048 .* Vz_m_day(:,ii,:) ./ days_per_month(ii);
end

In transient state the data are printed in multiple files. We refer to the file that contains the coordinates of the velocity points as the XYZ file. The following code puts all the data we need in a 2D array

nvel_p = size(MSH_ELEV,1)*size(MSH_ELEV,2);
XYZ_DATA = [...
    repmat(BCX_3310,4,1) ... % X
    repmat(BCY_3310,4,1) ... % X
    reshape(MSH_ELEV, nvel_p, 1) ... % Z
    zeros(nvel_p,1) ... % Proc ID
    repmat(DIAM,4,1) ... % DIAM 
    reshape(bsxfun(@rdivide, DIAM, av_thick), nvel_p, 1) ... % RATIO
];

We repeat a similar step for the velocities

VX_DATA = [];
VY_DATA = [];
VZ_DATA = [];
for ii = 1:size(Vx_m_day,3)
    VX_DATA = [VX_DATA; Vx_m_day(:,:,ii)];
    VY_DATA = [VY_DATA; Vy_m_day(:,:,ii)];
    VZ_DATA = [VZ_DATA; Vz_m_day(:,:,ii)];
end

Write the data of the last decade

• ASCII format

Write the XYZ file

prefix = fullfile(input_files_path, 'c2vsim_VelCLOUD_TR_');
fid = fopen([prefix 'XYZ_000.ich'],'w');
fprintf(fid, '%.3f %.3f %.3f %d %.2f %.2f\n', XYZ_DATA');
fclose(fid);

and the velocity files

mult = 1000000;
itimes = 385:504;
frmt = ['%.5f' repmat(' %.5f',1,length(itimes)-1) '\n'];

% VX
fid = fopen([prefix 'VX_000.ich'],'w');
fprintf(fid, frmt, mult*VX_DATA(:,itimes)');
fclose(fid);
% VY
fid = fopen([prefix 'VY_000.ich'],'w');
fprintf(fid, frmt, mult*VY_DATA(:,itimes)');
fclose(fid);
% VZ
fid = fopen([prefix 'VZ_000.ich'],'w');
fprintf(fid, frmt, mult*VZ_DATA(:,itimes)');
fclose(fid);

• HDF5 format For HDF5 format we write the data to a single file under different datasets as:

h5create([prefix '000.h5'],'/XYZDR',[size(XYZ_DATA,1) 5], 'Datatype','single');
h5write([prefix '000.h5'], '/XYZDR', XYZ_DATA(:,[1 2 3 5 6]));
h5create([prefix '000.h5'],'/PROC',[size(XYZ_DATA,1) 1], 'Datatype','uint32');
h5write([prefix '000.h5'], '/PROC', XYZ_DATA(:,4));
h5create([prefix '000.h5'],'/VX',[size(VX_DATA,1) size(VX_DATA,2)], 'Datatype','single');
h5write([prefix '000.h5'], '/VX', mult*VX_DATA);
h5create([prefix '000.h5'],'/VY',[size(VY_DATA,1) size(VY_DATA,2)], 'Datatype','single');
h5write([prefix '000.h5'], '/VY', mult*VY_DATA);
h5create([prefix '000.h5'],'/VZ',[size(VZ_DATA,1) size(VZ_DATA,2)], 'Datatype','single');
h5write([prefix '000.h5'], '/VZ', mult*VZ_DATA);

The graph file is the same in both steady and transient state.

The time step file for the last decade of C2VSim is:

Ndays = eomday(VelOut.YMD(itimes,1), VelOut.YMD(itimes,2));
TS = [0; cumsum(Ndays)];
fid = fopen(fullfile(input_files_path,'TimeStep10yrs.ich'),'w');
fprintf(fid,'%.1f\n', TS);
fclose(fid);

For this demo we have prepared two files. One for individual particles and for with wells.

In transient state we have to specify the release time. The above snippet sets the release time to 2161. In the TimeSteps variable we see that 2161 is between the 72 step which corresponds to SEP 2011. Therefore, these particles will be released the 30/SEP/2011.

In transient state the step size has to be at least as large as the time steps. The time steps are months therefore the steps will be generally very small. In fact, they are going to be so small that Euler method will give very accurate and smooth results.

The main configuration file for the transient state is the c2vsim_CLOUD_TR_config.ini and the velocity configuration is the c2vsim_CLOUD_TR_vf.ini.
Among the step configuration parameters the nStepsTime = 2 is important in transient state simulation. This parameter controls how many steps to take within a time step. Therefore since the difference between two steps is a month, setting this to 2 will force the code to take a step every 15 days.

In the following figure we compare the streamlines between steady (blue) and transient (red) state.

We can see that the steady state streamlines is a smoothed version of the transient state since the transient velocity field was repeated over and over.

The following plot shows the histograms of the streamline lengths and travel time. We can see that the length for the transient streamlines is are longer compared to steady, which is expected as their trajectories are more wiggly.
However the travel time is of the same order. This is because the steady state velocity is the average velocity. Although the transient streamlines are longer, the cumulative time of the low and high velocity segments is approximately equal to the cumulative average velocity.