SVSim cloud
In this example we will demonstrate how to prepare the input files for the SVSim model which is an IWFM based model for the nothern part of Central Valley.
Here we use the SVSim version 1.0 which can be fetched from https://data.cnra.ca.gov/dataset/svsim. The version used in this example is the June 20 2021. Before run the simulation we need to add the option to print the velocity field. Under the folder Simulation/Groundwater/ we will change a few parameters of the SVSim_Groundwater1973.dat which is the main groundwater input file for the model. First we will change the value of FACTVROU from 0.000022957 to 1000. Then we will change the value of UNITVROU from AC-FT/MON to 10^3FT/MON. The UNITVROU does not have any effect during hte simulation. It is only used to print the correct units in the output velocity file see also.
Next we set filenames fro the parameters VELOUTFL and VFLOWOUTFL
..\Results\SVSim_Velocity.out / VELOUTFL
..\Results\SVSim_VertVel.out / VFLOWOUTFL
After those changes we run SVSim with the new options. Before getting into the input files of Ichnos for the SVSim model, there are two output data from SVSim 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. Keep in mind that the velocity reader uses alot of memory and depending on the available memory it can be very slow.
svsim_path = fullfile('path','to','svsim');
Heads = readIWFM_headalloutput(fullfile(svsim_path,'Results','SVSim_GW_HeadAll.out'), 22642, 4, 505, 1);
VelOut = readIWFM_Velocity(fullfile(svsim_path,'Results','SVSim_Velocity.out'), 32537, 4, 1000);
Then we can save them.
save('SVSIM_GW_HeadAll','Heads')
save(fullfile('SVSIM_VELOUT'), 'VelOut','-v7.3');
For the domain outline we can read the elements and create an element shapefile which can then be disolved to a single polygon using a GIS software. For the SVSim the repository contains several gis layers. One of them is the SVSim_Boundary which contains the information we need.
Read the shapefile and split it into individual polygons.
sv_outline = shaperead(fullfile('path','to','SVSim_Boundary'));
[Xs, Ys] = polysplit(sv_outline.X, sv_outline.Y);
Next we print the domain file. Details on the format of the outline file can be found here
fid = fopen(fullfile('SV_outline.ich') ,'w');
for ii = 1:length(Xs)
[lat,lon] = projinv(projcrs(26910), Xs{ii,1}, Ys{ii,1});
[X3310, Y3310] = projfwd(projcrs(3310),lat, lon);
fprintf(fid, '%d %d\n',[length(X3310) ispolycw(X3310,Y3310)]);
fprintf(fid, '%.3f %.3f\n', [X3310' Y3310']');
end
fclose(fid);
In the above code we convert the coordinates from the EPSG 26910 which is the system that SVSim uses to 3310.
First we read the nodes from the SVSim file node file and convert them to 3310
ND = readIWFM_Nodes(fullfile(svsim_path,'PreProcessor','SVSim_Nodes.dat'));
[lat,lon] = projinv(projcrs(26910),ND(:,1), ND(:,2));
[ND_3310(:,1), ND_3310(:,2)] = projfwd(projcrs(3310),lat, lon);
Next we read the stratigraphy and calculate the bottom aquifer elevation in feet
STRAT = readIWFM_Stratigraphy(fullfile(svsim_path,'PreProcessor','SVSim_Stratigraphy.dat'),size(ND,1),9,108);
GSE = STRAT(:,2)*0.3048;
SVSimBot = GSE - sum(STRAT(:,3:end),2)*0.3048;
For the top we will use an average head of the last 5 years of the simulation. We load the head variable that we read previously.
load('path\to\SVSIM_GW_HeadAll.mat');
Next we calculate the average hydraulic head of the first layer for the last 5 years of the simulation (from OCT 2010 to SEP 2015)
AVHead = zeros(size(ND,1),1);
Ntimes = 0;
for ii = 446:505
AVHead = AVHead + Heads{ii,2}(:,1)*0.3038;
Ntimes = Ntimes + 1;
end
AVHead = AVHead./Ntimes;
Last we write the input file. :warning: It should be noted that we have to make sure that the elevation of bottom should be always less than the elevation of top.
fid = fopen(fullfile('SV_TopBot_CLOUD.ich'),'w');
fprintf(fid, 'CLOUD\n');
fprintf(fid, '%.2f %.2f\n', [6000, 3]);
fprintf(fid, '%.3f %.3f %.4f %.4f\n', [ND_3310 AVHead SVSimBot]');
fclose(fid);
In the second line the parameters correspond to radius and power of interpolation (see more). The radius should be set equal to a value so that when we draw a circle from any point of the domain with radius equal the value at least 3-5 points of the cloud points should be found.
One way to find such value is the following
nd_dists = zeros(size(ND_3310,1),6);
for ii = 1:size(ND_3310,1)
dst = sort(sqrt((ND_3310(ii,1) - ND_3310(:,1)).^2 + (ND_3310(ii,2) - ND_3310(:,2)).^2));
nd_dists(ii,:) = dst(1:6);
end
fprintf('%.2f ', max(nd_dists))
The last line displays the values
0.00 2796.34 3072.69 4937.90 5055.76 5979.34
If we select R = 3000 there will be areas in the domain where only 1 point can be found within the search circle. Here we set it equal to 6000 which means at least 5 nodes can be found within and distance R from any point of the domain.
Load the velocity data that we previously read
load('path\to\SVSIM_VELOUT.mat');
The saved variable is a structure that contains the coordinates of the 2D element barycenters where the velocity is defined. However, the coordinates are in ft and 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);
Next we calculate the elevation of the node barycenters. We will use the STRAT variable that contains the stratigraphy of the file. The stratigraphy file define elevation at the nodes.
LAYELEV_ND = zeros(size(ND_3310,1),10);
LAYELEV_ND(:,1) = STRAT(:,2)*0.3048;
for ii = 2:10
LAYELEV_ND(:,ii) = LAYELEV_ND(:,ii-1) - sum(STRAT(:,ii*2-1:ii*2)*0.3048,2);
end
Now we will loop through the elements and calculate the elevation of the barycenter, the dimateter and the ratio. In SVSim the node ID numbers do not coincide with the row numbers while node ids of the mesh files point to the node ids of the node file. To avoid this we renumber the mesh so that the node ids of the mesh elements point to the row number of the node ids. This will make the indexing much easier in the the next snippet.
MSH = readIWFM_Elements(fullfile(svsim_path,'PreProcessor','SVSim_Elements.dat'),23767,133);
MSH_renum = zeros(size(MSH,1),4);
for ii = 1:size(ND,1)
idx = find(MSH(:,2:5) == ND(ii,3));
MSH_renum(idx) = ii;
end
MSH_renum = [MSH(:,1) MSH_renum MSH(:,6)];
Now we can use the remumbered mesh to calculate the proprties mentioned previously
for ii = 1:size(MSH,1)
nv = 4;
if MSH_renum(ii,5) == 0
nv = 3;
end
for j = 1:9
t = mean(LAYELEV_ND(MSH_renum(ii,2:nv+1),j));
b = mean(LAYELEV_ND(MSH_renum(ii,2:nv+1),j+1));
MSH_ELEV(ii,j) = (t + b)/2;
av_thick(ii,j) = mean(LAYELEV_ND(MSH_renum(ii,2:nv+1),j) - LAYELEV_ND(MSH_renum(ii,2:nv+1),j+1));
end
DIAM(ii,1) = max(pdist([ND_3310(MSH_renum(ii,2:nv+1),1) ND_3310(MSH_renum(ii,2:nv+1),2)]));
end
Graph file or connectivity file shows how the nodes of the could are connected. It generally improves the interpolation of the velocity field.
First we create the 2D connection graph of the mesh. For each element we will make a list of the elements that each is connected to
msh_v = MSH_renum(:,2:5); % Isolate the mesh vertices
clear msh_graph
msh_graph{size(MSH_renum,1),1} = [];
for ii = 1:size(MSH_renum,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,9} = [];
for ii = 1:size(MSH_renum,1)
for jj = 1:9
if jj == 1
mesh3D_graph{ii,jj} = [jj*nel+ii; (jj-1)*nel+msh_graph{ii}; jj*nel+msh_graph{ii}];
elseif jj == 9
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, size(mesh3D_graph,1)*size(mesh3D_graph,2), 1);
Now we can print the graph file (see details on the file format).
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);
For the steady state simulation we will average the velocities of the last 5 years of the simulation (from OCT 2010 to SEP 2015)
for ii = 1:size(VelOut.VX,3)
Vx_av(:,ii) = mean(VelOut.VX(:,445:504,ii),2);
Vy_av(:,ii) = mean(VelOut.VY(:,445:504,ii),2);
Vz_av(:,ii) = mean(VelOut.VZ(:,445:504,ii),2);
end
Convert the velocity units to m/day. The original units of SVSim are FT/MONTH
Vx_av = Vx_av .* (0.3048/30.437);
Vy_av = Vy_av .* (0.3048/30.437);
Vz_av = Vz_av .* (0.3048/30.437);
Now we can reshape the data and print it to a file
nvel_p = size(MSH_ELEV,1)*size(MSH_ELEV,2);
VEL_DATA = [...
repmat(BCX_3310,9,1) ... % X
repmat(BCY_3310,9,1) ... % X
reshape(MSH_ELEV, nvel_p, 1) ... % Z
zeros(nvel_p,1) ... % Proc ID
repmat(DIAM,9,1) ... % DIAM
reshape(bsxfun(@rdivide, DIAM, av_thick), nvel_p, 1) ... % RATIO
reshape(Vx_av, nvel_p, 1) ... % VX
reshape(Vy_av, nvel_p, 1) ... % VY
reshape(Vz_av, nvel_p, 1) ... % VZ
];
Write the velocity 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);
or HDF5
fname = fullfile('svsim_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]));
Last we print the graph file we prepare above
fid = fopen(fullfile('svsim_VelCLOUD_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);
For the transient velocity field we will break the velocity components to different variables and print them under different files/datasets
First break into different variables and convert 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/datasets. 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,9,1) ... % X
repmat(BCY_3310,9,1) ... % X
reshape(MSH_ELEV, nvel_p, 1) ... % Z
zeros(nvel_p,1) ... % Proc ID
repmat(DIAM,9,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
Then we can write the files as ASCII XYZ file
prefix = fullfile('svsim_VelCLOUD_TR_');
fid = fopen([prefix 'XYZ_000.ich'],'w');
fprintf(fid, '%.3f %.3f %.3f %d %.2f %.2f\n', XYZ_DATA');
fclose(fid);
VX VY VZ files
itimes = 1:504;
frmt = ['%.10e' repmat(' %.10e',1,length(itimes)-1) '\n'];
% VX
fid = fopen([prefix 'VX_000.ich'],'w');
fprintf(fid, frmt, VX_DATA(:,itimes)');
fclose(fid);
% VY
fid = fopen([prefix 'VY_000.ich'],'w');
fprintf(fid, frmt, VY_DATA(:,itimes)');
fclose(fid);
% VZ
fid = fopen([prefix 'VZ_000.ich'],'w');
fprintf(fid, frmt, VZ_DATA(:,itimes)');
fclose(fid);
HDF5 format
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', VX_DATA);
h5create([prefix '000.h5'],'/VY',[size(VY_DATA,1) size(VY_DATA,2)], 'Datatype','single');
h5write([prefix '000.h5'], '/VY', VY_DATA);
h5create([prefix '000.h5'],'/VZ',[size(VZ_DATA,1) size(VZ_DATA,2)], 'Datatype','single');
h5write([prefix '000.h5'], '/VZ', VZ_DATA);
The time unit of the simulation is day as we converted the velocity to m/day. The velocity data in the previous snippet composed of 504 columns which is the number of timesteps. The time step file provides information on the number of days since the timestep 1. In the following snippet we use the itimes that was defined previously.
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);
The time step file contains one more element as it provides the start and ane time of each time step
T1---V1---T2---V2---T3---...---Tn---Vn---Tn+1
To demonstrate the tracing simulation we will generate 1000 particles distributed randomly inside the domain. To avoid the boundary effects the particles will be distributed within a buffer zone 5km from the outline. In the following snippet the first line creates a polygon object. The second line offsets the polygon inwords by 5000 and the 3rd line creates an interpolant for the top elevation. In the loop we generate 1000 points and their Z location is set 2 m below the top of the domain.
sv_outline.pp = polyshape(sv_outline.X, sv_outline.Y);
sv_buff = polybuffer(sv_outline.pp,-5000);
Ftop = scatteredInterpolant(ND(:,1),ND(:,2),AVHead);
particles = zeros(1000,6);
for ii = 1:1000
while 1
px = sv_outline.BoundingBox(1,1) + (sv_outline.BoundingBox(2,1) - sv_outline.BoundingBox(1,1))*rand;
py = sv_outline.BoundingBox(1,2) + (sv_outline.BoundingBox(2,2) - sv_outline.BoundingBox(1,2))*rand;
if sv_buff.isinterior(px, py)
break;
end
end
pz = Ftop(px,py)-2;
particles(ii,:) = [1 ii px py pz 0];
end
Before we print the file we have to convert the particle coordinates to EPSG:3310
[lat,lon] = projinv(projcrs(26910), particles(:,3), particles(:,4));
[particles(:,3), particles(:,4)] = projfwd(projcrs(3310),lat, lon);
fid = fopen('sv_particles.ich','w');
fprintf(fid,'%d %d %.3f %.3f %.3f %.1f\n', particles');
fclose(fid);
To run the code we prepare the main configuration file and the velocity file and then we run the model as
ichnos.exe -c svsim_cloud_config.ini
where svsim_cloud_config.ini is the main configuration file.
To conclude this demonstration we run a steady state simulation using RK45 and a transient simulation using Euler time step. Generally in transient groundwater simulations where the velocities are rather slow using a simple Euler scheme is capable to provide accurate results. In both simulations we set the MaxIterationsPerStreamline = 10000
We can read the results as
S_ss = readICHNOStraj(fullfile('output\sv_cloud_ss_iter_0000.h5'));
S_tr = readICHNOStraj(fullfile('output\sv_cloud_tr_iter_0000.h5'));
The following snippet plots the streamlines for the two simulations
cl_ord = colororder
figure()
clf
plot(X3310, Y3310,'LineWidth',2,'Color',cl_ord(1,:))
hold on
plot(particles(:,3),particles(:,4),'.','Color',cl_ord(2,:))
for ii = 1:size(S_ss)
if size(S_ss(ii,1).p,1) == 1
continue;
end
plot(S_ss(ii,1).p(:,1),S_ss(ii,1).p(:,2),'Color',cl_ord(3,:))
end
for ii = 1:size(S_tr)
if size(S_tr(ii,1).p,1) == 1
continue;
end
plot(S_tr(ii,1).p(:,1),S_tr(ii,1).p(:,2),'Color',cl_ord(4,:))
end
axis equal off
The left figure shows all the streamlines and the right zooms into one area. We can see that the yellow streamlines (steady state) are smoother and complete. They are smoother because they are based on a steady state flow field and complete because we use RK45 which allows to trace the full paths. The path trajectories for the transient state simulation appear jagged because the velocity field exhibit monthly changes and they are not complete because the 10000 limit didn't allow to delinieate full paths.