Embedded Files
Introduction
The mathematical model we developed weds numerical analysis and computer science as applied practically to measure and predict the dynamically changing topography, salinity, water level, and exposed shoreline area of the Salton Sea over 50 years starting from 2003.
Background
In the early 20th century the Imperial Canal, which was developed to supply water to the Imperial Valley from the Colorado River, was breached and for approximately a year and a half, the Colorado River flowed freely in to the Salton sink.
This newly formed body of water is currently the largest man-made inland body of water in the state of California and is commonly called the Salton Sea. As of 2005, the Salton Sea is 15 miles by 35 miles with a shoreline elevation of approximately 228.7 feet below sea level. The main sources of inflows to the Salton Sea are the following: Alamo River, New River, Whitewater River, and direct agricultural drainage from the Imperial and Coachella valleys
In 2003, the Imperial Irrigation District (IID), which has rights to a large portion of the Colorado River, has agreed to provide more water to the San Diego and Los Angeles areas under the California Quantification Settlement Agreement (QSA). Starting at the end of 2017 this agreement will reduces inflows to from the Colorado River to supply more water to the San Diego and Los Angeles area. This reduction will not only affect water levels and salinity, but impact the environment, avian, and aquatic wildlife.
The impact of reducing inflow to the Salton Sea will have devastating effects on local habitat. The Salton Sea is located along one of the most important flyways in North America, providing critical habitat for more than 400 species of resident and migrating birds. As water elevation in the Salton Sea continues to decline and salinity continues to increase, the result is a decline in the habitat the birds and fish depend on and an eventual collapse of the fisheries. With no outflows aside from evaporation, the salinity in the Salton Sea is reported to be 33 % greater than that of the ocean.
Design Procedures
Entering the lake into MatLab®
We first began by constructing a three dimensional model of the of the Salton sea using a CAD program, SketchUp® by importing a two dimensional contour map of the Salton Sea then keeping 1 km in the contour map equal to 1 m in SketchUp®. This is important otherwise the model would be difficult to work with because the length and width of the lake are very large compared its depth. Using the line-tool we traced the contours (see Figure 2).
The contours were extruded to their appropriate depths in meters in SketchUp® and each of the contours lofted to their closest contour line, Figure 3. Once the general shape of the lake was made we needed have this same shape made into a grid with equally spaced x and y varying z coordinates (depth). A sketchUp® add-on was used to create this. The resulting lake is shown in Figure 4. This is important because the resulting intersections of the grid will act as the finite elements in the MatLab® simulation. The lake’s x and y coordinates were then scaled 1000 times so lake was the correct scale with each grid spacing being 1 km apart.
Finally, the lake was exported as an .stl file from SketchUp®. .stl is a universal CAD export format that simply exports the coordinate locations and face locations of each hidden line intersection and face of the model respectively from the created model. Once this model was exported the file was placed into the MatLab file folder along with the stlread MatLab® function. The stlread function imports .stl files into MatLab® and places all the hidden line intersections into a (n,3) matrix with n being six times the number of intersections. This is because the .stl file stores six values for each intersection this will be important to correct when writing the simulations or else the area and volumes will be six times larger than the real value. The 1st column is the x coordinate, 2nd column the y value and the 3rd column the z value (or depth). The model showing the lake at its initial conditions entered into MatLab® is shown in Figure 5. Each one of these points will represent our finite element with an area of 1 km2 or 1*106 m2.
Now, using this model in Matlab® we can perform calculations of changes in area, volume and salinity for each element in the model. Using this method for entering the lake coordinates we received less than a 5 % error in our initial volume (5.6*106 ac-ft) and less than 0.5 % error in our initial area (894 km2) compared the known volume and area of the lake, 6*106 ac-ft and 889 km2 respectivly.6
Assumptions
Before we began to calculate the changes in the lake by cuts in inflows we made the following assumptions:
Water surface area would change slowly enough to allow for recalculation of the lake surface on a yearly basis.
We assumed that the volume within the interpolated contour “mesh” generated in SketchUp® (See Figure 4) would be adequately modeled with Riemann’s sums of height equal to each z-data point on the mesh.
Our analysis assumes that average sun exposure over a year period adequately models the driving force of the evaporation
Since rain directly entering the lake is minimal compared to evaporation we will neglect it and assume most of the rain enters in from the inlet streams.
We assumed no salt left the Salton Sea basin.
The average salt content, inflows and outflows can be used to represent the instantaneous salt content, inflows and outflows and assume they do not vary significantly on a yearly basis.
*more assumptions that were made about the calculations and equations used are mentioned thought this paper.
We the applied these assumptions to calculate the changes in the lake on a yearly basis.
Equations and calculation methods
Once the lakes x, y and z coordinates have been entered into MatLab® iteration were made on a yearly basis to calculate the desired results.
The evaporation rate was assumed constant each year and was calculated using Hamon’s equation
Where E = evaporation of day t (mm day-1), Ht = average number of daylight hours per day during the month in which day t falls, es = saturated vapor pressure at temperature T (kPa) (see Penman’s equation below), Tt = temperature, day t (ºC). The temperatures we used were the daily mean temperatures for each month based on historical data7. Ht can be calculated by using the maximum number of daylight hours on day t, Nt. that is
Penman’s equation
These values were plotted in excel and the value of E was summed up for each day to obtain the overall E. The overall E obtained was 1.143 m year-1. We assumed this number multiplied by the surface area of the lake will give the volume of water evaporated each year. Also, the addition of the inflows multiplied by one year will give the volume entering the lake each year. Thus a new volume can be calculated for each year by the following equation
This equation was implemented into a while loop in MatLab® for each year iteration, i. If statements were used to implement the inflows and cut offs for each year. All the z elements greater than zero for that year are summed and subtracted them from the volume in the next iteration. This is used as a correction of elements whom have passed the zero depth point. Next all values greater than zero are set to zero and kept that way to prevent them to be used in later iterations.
The total salt in the lake can be calculated in a similar method. The salt coming in is calculated from the rate of water inflow multiplied by the salt concentration. The concentration of salt in the sea was calculated in each iteration by:
With Salt(out) assumed to be zero. Vi is the volume at that iteration. The exposed surface area of the lake was calculated by the number of z values or depth above zero. We did this by making the initial shore line elevation equal to zero in our simulation. We used -0.00000001 instead of zero to account for .stl to MatLab® rounding errors.
This equation of area was used to calculate the amount of water leaving though the ground and evaporation. We assumed that a percent change in area was equal to the same percent change in evaporation rate and percentage change of water leaving though the ground.
The final depth of the lake was found by finding the minimum z value in the lake coordinate matrix for each iteration. The change in elevation was found by knowing the initial shore line elevation (-235ft) and subtracting the change in lake depth.
This same method can be applied to all three scenarios given for the lake. In scenario 2.1 subtracting the inflow reductions from Mexicali and the accompanying salt inflow. In scenario 2.2 an addition cut off was made from the power plant in Mexicali and the accompanying salt was also removed from being added to the Salton Sea. Since the salt concentration given to the power plant was not given we assumed it to be the same concentration of salt entering the US-Mexico border, 3.5 tons/ac-ft. order of operations and code in the MatLab® is crucial to obtain accurate lake change calculations. The results are depicted in animations created in the MatLab® program and can be called by entering LakeSim1, LakeSim21 and LakeSim22 in the command line, representing scenario 1, scenario 2.1 and scenario 2.2 respectively . These must be in presence of the Lake.stl and stlread function in the same folder.
Results & Discussion
Using MATLAB we generated a three dimensional model of the contours of the topography of the Salton Sea. The initial water elevation, exposed shoreline and salt concentration was determined to be 13.0620 m, 235 ft. below sea level, 54 g L-1, respectively.
The results of our analysis for scenario 1 shown in Figure 6, paint a bleak future for the Salton Sea. According to our model, in fifty years the volume of the lake will be 61 % of the original volume from 2003. Likewise, by 2053, the surface area of the lake will be 82 % of the original volume and the salinity increased by approximately a factor of 2.22 to a final salinity of 120 g L-1. The exposed shore area increases to approximately 150 Km2.
The results of our analysis for scenario 2.1 shown in Figure 7, depict that in fifty years the volume of the lake will be 57.16 % of the original volume from 2003. Likewise, by 2053, the surface area of the lake will be 81.31 % of the original volume and the salinity will have increased by approximately 10 g L-1 from scenario 1 to 130 g L-1. The exposed shore area increases to approximately 168 Km2.
The results of our analysis for scenario 2.2 shown in Figure 8, depicts in fifty years the volume of the lake will be 54.04 % of the original volume from 2003. Likewise, by 2053, the surface area of the lake will be 79.86% of the original volume and the salinity will have increased by approximately 15 g L-1 from scenario 1 to 135 g L-1. The exposed shore area increases to approximately 180 Km2.
Alternate Solution:
Over the years, numerous solutions have been proposed to save the Salton Sea.
The Salton Sea Authority was created by the state of California "for the purpose of ensuring the beneficial uses of the Salton Sea." As such, the SSA is uniquely qualified to recommend courses of action regarding the lake and its ecosystem.
In addition to its official recommendations, the Salton Sea Authority provided a report in May of 2007 on alternate solutions to the maintenance of the Salton Sea ecosystem. In this report, the SSA outlined 8 potential alternate plans to its Salton Sea Ecosystem Restoration Program. The fifth alternative reported, which the SSA designated its “preferred” alternative, includes a Saline Habitat Complex in the northern and southern sea bed, a Marine Sea that extends around the northern shoreline, and Air Quality Management facilities to reduce particulate emissions. These features are shown in Figure 9.
In analyzing the alternate recommendations of the SSA, we must come back to our modeling and analysis. Regardless of the benefits of brine sinks and saline habitat complexes, and even air quality management programs, the simple fact remains that lakes need water. Our analysis is unequivocal in demonstrating that the Salton Sea will dry up if inflows are reduced. No amont of dressing the sea up will change that.
Conclusion & Recommendations
If the QSA is implemented, 30 million acre-ft of water over a period of 75 years will be diverted from farms in the Imperial Valley to urban areas in Southern California. Beginning in 2018, this will have a major impact on the Salton Sea and the surrounding ecosystem. Agricultural runoff that sustains the lake will decrease, and as our model predicts, the Salton Sea will dry up.
Having analyzed the result of reducing inflows to the lake, we must now address the ethical question of what to do with this knowledge.
We know that the Salton Sea is currently the largest lake in California. A lake situated in a desert environment, it uniquely provides for life in the region. Furthermore, the Salton Sea an important link in the bird migratory path known as the Pacific Flyway.
We also know that the entire state of California is experiencing a severe drought. We cannot sacrifice the lives of humans for those of birds. Therefore, our ethical responsibility as scientists is to report the impending death of the Salton Sea, but recommend that the QSA be acted upon.
References
Herting, Alex, Tim Farmer, and Jordan Evans. "Mapping of the Evaporative Loss from Elephant Butte Reservoir Using Remote Sensing and GIS Technology." NMSU CAGE (2004): n. pag. New Mexico State University. Web. 3 Dec. 2014. <http://wrri.nmsu.edu/research/rfp/studentgrants03/reports/herting.pdf>.
"Quantification Settlement Agreement." San Diego County Water Authority. N.p., n.d. Web. 10 Dec.2014. <http://www.sdcwa.org/quantification-settlement-agreement>.
"Salton Sea Ecosystem Restoration Program." Department of Water Resources. N.p., n.d. Web. 10 Dec. 2014. <http://www.water.ca.gov/saltonsea/>.
"Salton Sea 101." State of California. California Department of Parks and Recreation, n.d. Web. 10 Dec. 2014. <http://www.parks.ca.gov/?page_id=21264>.
Salton Sea Authority. N.p., n.d. Web. 10 Dec. 2014. <http://www.saltonsea.ca.gov>."WATER BOND PRIORITIES: California’s responsibility and greatest opportunity to revitalize a dying ecosystem." Salton Sea Authority. N.p., n.d. Web. 10 Dec. 2014. <http://saltonsea.ca.gov/ Portals/0/Future/Salton%20Sea%20Funding%20and%20Feasibility%20Review%20Work%20Plan.pdf>.
“Salton Sea Geography”, Salton Sea Authority, Web 10/17/2014 http://www.saltonsea.ca.gov/saltonsea/geography/
MSN Weather, Windows 8 weather App, Historical Weather, Brawley, California, referenced from “http://www.wdtinc.com/windows-8-consumer.com”
Open Drop Down to See Code
Open Drop Down to See Code
Appendix/MatLab® Code
.STL File Reader Function Code sltread();
function varargout = stlread(file)
% STLREAD imports geometry from an STL file into MATLAB.
% FV = STLREAD(FILENAME) imports triangular faces from the ASCII or binary
% STL file idicated by FILENAME, and returns the patch struct FV, with fields
% 'faces' and 'vertices'.
%
% [F,V] = STLREAD(FILENAME) returns the faces F and vertices V separately.
%
% [F,V,N] = STLREAD(FILENAME) also returns the face normal vectors.
%
% The faces and vertices are arranged in the format used by the PATCH plot
% object.
% Copyright 2011 The MathWorks, Inc.
if ~exist(file,'file')
error(['File ''%s'' not found. If the file is not on MATLAB''s path' ...
', be sure to specify the full path to the file.'], file);
end
fid = fopen(file,'r');
if ~isempty(ferror(fid))
error(lasterror); %#ok
end
M = fread(fid,inf,'uint8=>uint8');
fclose(fid);
[f,v,n] = stlbinary(M);
%if( isbinary(M) ) % This may not be a reliable test
% [f,v,n] = stlbinary(M);
%else
% [f,v,n] = stlascii(M);
%end
varargout = cell(1,nargout);
switch nargout
case 2
varargout{1} = f;
varargout{2} = v;
case 3
varargout{1} = f;
varargout{2} = v;
varargout{3} = n;
otherwise
varargout{1} = struct('faces',f,'vertices',v);
end
end
function [F,V,N] = stlbinary(M)
F = [];
V = [];
N = [];
if length(M) < 84
error('MATLAB:stlread:incorrectFormat', ...
'Incomplete header information in binary STL file.');
end
% Bytes 81-84 are an unsigned 32-bit integer specifying the number of faces
% that follow.
numFaces = typecast(M(81:84),'uint32');
%numFaces = double(numFaces);
if numFaces == 0
warning('MATLAB:stlread:nodata','No data in STL file.');
return
end
T = M(85:end);
F = NaN(numFaces,3);
V = NaN(3*numFaces,3);
N = NaN(numFaces,3);
numRead = 0;
while numRead < numFaces
% Each facet is 50 bytes
% - Three single precision values specifying the face normal vector
% - Three single precision values specifying the first vertex (XYZ)
% - Three single precision values specifying the second vertex (XYZ)
% - Three single precision values specifying the third vertex (XYZ)
% - Two unused bytes
i1 = 50 * numRead + 1;
i2 = i1 + 50 - 1;
facet = T(i1:i2)';
n = typecast(facet(1:12),'single');
v1 = typecast(facet(13:24),'single');
v2 = typecast(facet(25:36),'single');
v3 = typecast(facet(37:48),'single');
n = double(n);
v = double([v1; v2; v3]);
% Figure out where to fit these new vertices, and the face, in the
% larger F and V collections.
fInd = numRead + 1;
vInd1 = 3 * (fInd - 1) + 1;
vInd2 = vInd1 + 3 - 1;
V(vInd1:vInd2,:) = v;
F(fInd,:) = vInd1:vInd2;
N(fInd,:) = n;
numRead = numRead + 1;
end
end
function [F,V,N] = stlascii(M)
warning('MATLAB:stlread:ascii','ASCII STL files currently not supported.');
F = [];
V = [];
N = [];
end
% TODO: Change the testing criteria! Some binary STL files still begin with
% 'solid'.
function tf = isbinary(A)
% ISBINARY uses the first line of an STL file to identify its format.
if isempty(A) || length(A) < 5
error('MATLAB:stlread:incorrectFormat', ...
'File does not appear to be an ASCII or binary STL file.');
end
if strcmpi('solid',char(A(1:5)'))
tf = false; % ASCII
else
tf = true; % Binary
end
end
START LAKE SIMULATION SCENARIO 1 LakeSim1
i=0;
%these set inital conditions for Simulation
ShoreEO=-71.628; %inital shore line elevation
SEM=zeros(50,2); %Shore elevation matrix
numDay=1;
input=0;
evap=1.139;
out=1.0;
y=0; %year
Salt_in=1.06933*10^11;
Salt_inO=1.06933*10^11;
Salt_out=0;
Salt_Conc=0;
Salt=0;
Salinity=0;
SM=zeros(50,2); %50 is the number of years used to make plot of salinity
SA=zeros(50,2);
m=1;
n=1;
r=1;
ground_out=0;
p=1;
q=1;
d=2003;
Vcor=0;
S=0;
LakeCord=0;
Volume=0;
LakeAreaO=0;
acer_feet=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
LakeCord=stlread('Lake.stl'); % Creates Matrix in MatLab from .stl File
LakeCord=LakeCord.vertices; % Makes vertices in .stl into XYZ matrix
volume= -(1000000/6)*sum(LakeCord(:,3)); % Calculates the Inital Volume of the Lake
volumeO=volume;
Lake_DepthO=-min(LakeCord(:,3));
Lake_Depth= -min(LakeCord(:,3));
S=size(LakeCord); % Finds size of matrix for loop calculations
a=S(1,1)+1; % 'a' designates number of loops
while r<a % Calculates the Original Lake Area by number of vertices below zero
if(LakeCord(r,3)< -0.00000000001)
LakeAreaO=LakeAreaO+1;
end
r = r + 1;
end
r=1;
LakeAreaO=LakeAreaO*1000000/6;
LakeArea=LakeAreaO;
ground_outO=1.43697*10^9-evap*LakeAreaO;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Salt=volume*54;
Salt_out=ground_outO*Salt/volume;
Salinity=Salt/volume;
%Calculates the initial amount of salt in the lake 54kg/m^3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Start Simulation %%%%%%%%%%%%%%%%%%%%%%%
while i < 50
if y<=10
input=(1.43697*10^9);
input=input-y*(1.23348*10^7);
end
if (10<y)&&(y<16)
input=(1.31362*10^9);
end
if (16<=y)&&(y<17)
input=(1.27662*10^9);
end
if (17<=y)&&(y<18)
input=(1.23961*10^9);
end
if (18<=y)&&(y<19)
input=(1.20261*10^9);
end
if y>=19
input=(1.19027*10^9);
end
%%%%%%%%%%%%%%%%%%%%%%%%%Total Reductions%%%%%%%%%%%%%%%%%%%%
Salt_in = input*2.2063;
y=y+1;
%%%%%%%%%%%%%%%%%%%% Plotting %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
set(gcf,'Position',[100 0 1500 900])
if i<50
SM(i+1,1)=Salinity;
SM(i+1,2)=i;
end
if i<50
SEM(i+1,1)=(ShoreEO-(Lake_DepthO-Lake_Depth))*3.28084;
SEM(i+1,2)=i;
end
if i<50
SA(i+1,1)=(LakeAreaO-LakeArea)/1000000;
SA(i+1,2)=i;
end
subplot(1,3,1)
[ax,p1,p2]= plotyy(SM(:,2),SM(:,1),SEM(:,2),SEM(:,1),'plot', 'plot');
set(ax(1),'Ylim',[50 150])
set(ax(2),'Ylim',[-250 -230])
set(ax(1),'xlim',[0 50]);
set(ax(2),'xlim',[0 50]);
set(ax(1),'YTick',50:10:150)
set(ax(2),'YTick',-250:1:-230)
set(p1,'marker','.')
set(p2,'marker','.')
grid(ax(1),'on')
ylabel('Salinity (g/L)')
set(get(ax(2),'Ylabel'),'string','Shore Line Elevation (ft below sea level)')
subplot(1,3,2);
plot3(LakeCord(:,1),LakeCord(:,2),LakeCord(:,3),'.');
axis([0 50000 0 50000 -13.5 0]);
xlabel('meters')
ylabel('meters')
zlabel('depth of lake (meters)')
title(d);
d=d+1;
pause(0.005);
subplot(1,3,3);
plot(SA(:,2),SA(:,1),'.');
grid('on')
title('Exposed Shore Area')
xlabel('year')
ylabel('square kilometers (km^2)')
axis([0 50 0 200])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Lake Change Calculations %%%%%%%%%%%%%%%%%%%%%%%%%%%%
out=0;
LakeArea=0;
while m<a % Calculates the Lake water Area by number of vertices below zero
if(LakeCord(m,3)< -0.0000000001)
LakeArea=LakeArea+1;
end
m = m + 1;
end
m=1;
LakeArea=LakeArea*(1000000/6);
ground_out=ground_outO*(LakeArea/LakeAreaO);
out=evap*LakeArea+ground_out;
Salt=Salt+Salt_in;
volume = -(1000000/6)*sum(LakeCord(:,3));
while p < a % Calculates drop in water height of each vertical element
if (LakeCord(p,3) <= -0.0000000001)
LakeCord(p,3) = LakeCord(p,3) - ((input-out)/LakeArea)-Vcor/LakeArea;
end
Vcor=0;
p= p + 1;
end
p=1;
q=1;
while q < a % Corrects the volume if the change in volume falls below an element.
if LakeCord(q,3) > 0
Vcor = Vcor + (1000000/6)*LakeCord(q,3);
end
q=q+1;
end
subLakeCord=LakeCord(:,3); % Assigns 0 to all non negative depth values to the lake to calculate new area
subLakeCord(subLakeCord>0)=0;
LakeCord(:,3)=subLakeCord;
LakeCord;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Lake_Depth= -min(LakeCord(:,3));
Percent_of_Original_Volume = (volume/volumeO)*100;
acer_feet=0.00081071*volume;
LakeArea;
y;
i=i+1;
input;
out;
Salinity=Salt/volume;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Program Outputs%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Oringinal_Area_of_Lake_square_kilometers=(LakeAreaO/1000000)
exposed_Shoreline_Square_kilometers =(LakeAreaO-LakeArea)/1000000% Original Area of lake
Original_Acer_feet_of_Lake=volumeO/1233.48184 %Original Volume of lake (2006)
Final_Acer_Feet=acer_feet
volume; % volume of Lake
LakeArea; % Exposed surface area of lake
Original_Lake_Depth_meters=Lake_DepthO
Lake_Depth_meters=Lake_Depth % Lake Depth and the end of the Simulation
Percent_of_Original_Volume
acer_feet;
START LAKE SIMULATION SCENARIO 2.1 LakeSim21
%these set inital conditions for Simulation
SA=zeros(50,2);
mexicali=0;
powerplant = 0;
ShoreEO=-71.628; %inital shore line elevation
SEM=zeros(50,2); %Shore elevation matrix
numDay=1;
input=0;
evap=1.139;
i=0;
out=1.0;
y=0; %year
Salt_in=1.06933*10^11;
Salt_inO=1.06933*10^11;
Salt_out=0;
Salt_Conc=0;
Salt=0;
Salinity=0;
SM=zeros(50,2); %50 is the number of years used to make plot of salinity
m=1;
n=1;
r=1;
ground_out=0;
p=1;
q=1;
d=2003;
Vcor=0;
S=0;
LakeCord=0;
Volume=0;
LakeAreaO=0;
acer_feet=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
LakeCord=stlread('Lake.stl'); % Creates Matrix in MatLab from .stl File
LakeCord=LakeCord.vertices; % Makes vertices in .stl into XYZ matrix
volume= -(1000000/6)*sum(LakeCord(:,3)) % Calculates the Inital Volume of the Lake
volumeO=volume;
Lake_DepthO=-min(LakeCord(:,3));
Lake_Depth=-min(LakeCord(:,3));
%%%%%%%%%%%%% Finds size of matrix for loop calculations %%%%%%%%%%%%%%%%%
S=size(LakeCord);
a=S(1,1)+1; % 'a' designates number of loops
%%%%% Calculates the Original Lake Area by number of vertices below zero%%%
while r<a
if(LakeCord(r,3)< -0.00000000001)
LakeAreaO=LakeAreaO+1;
end
r = r + 1;
end
r=1;
LakeAreaO=LakeAreaO*1000000/6;
LakeArea=LakeAreaO;
ground_outO=1.43697*10^9-evap*LakeAreaO;
%%%%%%%%% Calculates the initial amount of salt in the lake 54kg/m^3 %%%%
Salt=volume*54;
Salinity=Salt/volume;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Start Simulation %%%%%%%%%%%%%%%%%%%%%%%
while i < 50 % Stop scrpit once equilibrium is acheived
if y<=10
input=(1.43697*10^9);
input=input-y*(1.23348*10^7);
end
if (10<y)&&(y<16)
input=(1.31362*10^9);
end
if (16<=y)&&(y<17)
input=(1.27662*10^9);
end
if (17<=y)&&(y<18)
input=(1.23961*10^9);
end
if (18<=y)&&(y<19)
input=(1.20261*10^9);
end
if y>=19
input=(1.19027*10^9);
end
%%%%%%%%%%%%%%%%%%%%%% Mexicali Reductions %%%%%%%%%%%%%%%%%%%%%%%%%
if y>=3
mexicali=(1.85022*10^-7);% 1.85022*10^7m^3/year = 15,000 ac-ft/year reduction from mexicali
end
if y>=11
mexicali=(2.7762*10^7);% 2.7762*10^7m^3/year = 22,507 ac-ft/year reduction from mexicali
end
%%%%%%%%%%%%%%%%%%%%%%% PowerPlant %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%if y>=12
% powerplant=(1.97375*10^7);
%end
%%%%%%%%%%%%%%%%%%%%% Total Reductions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Salt_in = (input*2.2063) - (mexicali*0.9665);%- (powerplant*2.2063);
input=input-mexicali;%-powerplant;
y=y+1;
%%%%%%%%%%%%%%%%%%%%%%%% Plotting %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
set(gcf,'Position',[100 0 1500 900])
if i<50
SM(i+1,1)=Salinity;
SM(i+1,2)=i;
end
if i<50
SEM(i+1,1)=(ShoreEO-(Lake_DepthO-Lake_Depth))*3.28084;
SEM(i+1,2)=i;
end
if i<50
SA(i+1,1)=(LakeAreaO-LakeArea)/1000000;
SA(i+1,2)=i;
end
subplot(1,3,1)
[ax,p1,p2]= plotyy(SM(:,2),SM(:,1),SEM(:,2),SEM(:,1),'plot', 'plot');
set(ax(1),'Ylim',[50 150])
set(ax(2),'Ylim',[-250 -230])
set(ax(1),'xlim',[0 50]);
set(ax(2),'xlim',[0 50]);
set(ax(1),'YTick',50:10:150)
set(ax(2),'YTick',-250:1:-230)
set(p1,'marker','.')
set(p2,'marker','.')
grid(ax(1),'on')
ylabel('Salinity (g/L)')
set(get(ax(2),'Ylabel'),'string','Shore Line Elevation (ft below sea level)')
subplot(1,3,2);
plot3(LakeCord(:,1),LakeCord(:,2),LakeCord(:,3),'.');
axis([0 50000 0 50000 -13.5 0]);
xlabel('meters')
ylabel('meters')
zlabel('depth of lake (meters)')
title(d);
d=d+1;
pause(0.005);
subplot(1,3,3);
plot(SA(:,2),SA(:,1),'.');
grid('on')
title('Exposed Shore Area')
xlabel('year')
ylabel('square kilometers (km^2)')
axis([0 50 0 200])
%%%%%%%%%%%%%%%%%%%%%%%%%% Lake Calculations %%%%%%%%%%%%%%%%%%%%%%%%%%%%
out=0;
LakeArea=0;
while m<a % Calculates the Lake water Area by number of vertices below zero
if(LakeCord(m,3)< -0.0000000001)
LakeArea=LakeArea+1;
end
m = m + 1;
end
m=1;
LakeArea=LakeArea*(1000000/6);
ground_out=ground_outO*(LakeArea/LakeAreaO);
out=evap*LakeArea+ground_out;
Salt=Salt+Salt_in;
volume = -(1000000/6)*sum(LakeCord(:,3));
while p < a % Calculates drop in water height of each vertical element
if (LakeCord(p,3) <= -0.0000000001)
LakeCord(p,3) = LakeCord(p,3) - ((input-out)/LakeArea)-Vcor/LakeArea;
end
Vcor=0;
p= p + 1;
end
p=1;
q=1;
while q < a % Corrects the volume if the change in volume falls below an element.
if LakeCord(q,3) > 0
Vcor = Vcor + (1000000/6)*LakeCord(q,3);
end
q=q+1;
end
subLakeCord=LakeCord(:,3); % Assigns 0 to all non negative depth values to the lake to calculate new area
subLakeCord(subLakeCord>0)=0;
LakeCord(:,3)=subLakeCord;
LakeCord;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Lake_Depth= -min(LakeCord(:,3));
Percent_of_Original_Volume = (volume/volumeO)*100;
acer_feet=0.00081071*volume;
LakeArea;
y;
i=i+1;
input;
out;
Salinity=Salt/volume;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Program Outputs%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Oringinal_Area_of_Lake_square_kilometers=(LakeAreaO/1000000)
exposed_Shoreline_Square_kilometers =(LakeAreaO-LakeArea)/1000000% Original Area of lake
Original_Acer_feet_of_Lake=volumeO/1233.48184 %Original Volume of lake (2006)
Final_Acer_Feet=acer_feet
volume; % volume of Lake
Final_Lake_Area_square_kilometers=LakeArea/1000000 % Exposed surface area of lake
Original_Lake_Depth_meters=Lake_DepthO
Lake_Depth_meters=Lake_Depth % Lake Depth and the end of the Simulation
Percent_of_Original_Volume
acer_feet;
START LAKE SIMULATION SCENARIO 2.2 LakeSim22
%these set inital conditions for Simulation
SA=zeros(50,2);
mexicali=0;
powerplant = 0;
ShoreEO=-71.628; %inital shore line elevation
SEM=zeros(50,2); %Shore elevation matrix
numDay=1;
input=0;
evap=1.139;
i=0;
out=1.0;
y=0; %year
Salt_in=1.06933*10^11;
Salt_inO=1.06933*10^11;
Salt_out=0;
Salt_Conc=0;
Salt=0;
Salinity=0;
SM=zeros(50,2); %50 is the number of years used to make plot of salinity
m=1;
n=1;
r=1;
ground_out=0;
p=1;
q=1;
d=2003;
Vcor=0;
S=0;
LakeCord=0;
Volume=0;
LakeAreaO=0;
acer_feet=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
LakeCord=stlread('Lake.stl'); % Creates Matrix in MatLab from .stl File
LakeCord=LakeCord.vertices; % Makes vertices in .stl into XYZ matrix
volume= -(1000000/6)*sum(LakeCord(:,3)) % Calculates the Inital Volume of the Lake
volumeO=volume;
Lake_DepthO=-min(LakeCord(:,3));
Lake_Depth=-min(LakeCord(:,3));
%%%%%%%%%%%%% Finds size of matrix for loop calculations %%%%%%%%%%%%%%%%%
S=size(LakeCord);
a=S(1,1)+1; % 'a' designates number of loops
%%%%% Calculates the Original Lake Area by number of vertices below zero%%%
while r<a
if(LakeCord(r,3)< -0.00000000001)
LakeAreaO=LakeAreaO+1;
end
r = r + 1;
end
r=1;
LakeAreaO=LakeAreaO*1000000/6;
LakeArea=LakeAreaO;
ground_outO=1.43697*10^9-evap*LakeAreaO;
%%%%%%%%% Calculates the initial amount of salt in the lake 54kg/m^3 %%%%
Salt=volume*54;
Salinity=Salt/volume;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Start Simulation %%%%%%%%%%%%%%%%%%%%%%%
while i < 50 % Stop scrpit once equilibrium is acheived
if y<=10
input=(1.43697*10^9);
input=input-y*(1.23348*10^7);
end
if (10<y)&&(y<16)
input=(1.31362*10^9);
end
if (16<=y)&&(y<17)
input=(1.27662*10^9);
end
if (17<=y)&&(y<18)
input=(1.23961*10^9);
end
if (18<=y)&&(y<19)
input=(1.20261*10^9);
end
if y>=19
input=(1.19027*10^9);
end
%%%%%%%%%%%%%%%%%%%%%% Mexicali Reductions %%%%%%%%%%%%%%%%%%%%%%%%%
if y>=3
mexicali=(1.85022*10^-7);% 1.85022*10^7m^3/year = 15,000 ac-ft/year reduction from mexicali
end
if y>=11
mexicali=(2.7762*10^7);% 2.7762*10^7m^3/year = 22,507 ac-ft/year reduction from mexicali
end
%%%%%%%%%%%%%%%%%%%%%%% PowerPlant %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if y>=12
powerplant=(1.97375*10^7);
end
%%%%%%%%%%%%%%%%%%%%% Total Reductions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Salt_in = (input*2.2063) - (mexicali*0.9665)- (powerplant*2.2063);
input=input-mexicali-powerplant;
y=y+1;
%%%%%%%%%%%%%%%%%%%%%%%% Plotting %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
set(gcf,'Position',[100 0 1500 900])
if i<50
SM(i+1,1)=Salinity;
SM(i+1,2)=i;
end
if i<50
SEM(i+1,1)=(ShoreEO-(Lake_DepthO-Lake_Depth))*3.28084;
SEM(i+1,2)=i;
end
if i<50
SA(i+1,1)=(LakeAreaO-LakeArea)/1000000;
SA(i+1,2)=i;
end
subplot(1,3,1)
[ax,p1,p2]= plotyy(SM(:,2),SM(:,1),SEM(:,2),SEM(:,1),'plot', 'plot');
set(ax(1),'Ylim',[50 150])
set(ax(2),'Ylim',[-250 -230])
set(ax(1),'xlim',[0 50]);
set(ax(2),'xlim',[0 50]);
set(ax(1),'YTick',50:10:150)
set(ax(2),'YTick',-250:1:-230)
set(p1,'marker','.')
set(p2,'marker','.')
grid(ax(1),'on')
ylabel('Salinity (g/L)')
set(get(ax(2),'Ylabel'),'string','Shore Line Elevation (ft below sea level)')
subplot(1,3,2);
plot3(LakeCord(:,1),LakeCord(:,2),LakeCord(:,3),'.');
axis([0 50000 0 50000 -13.5 0]);
xlabel('meters')
ylabel('meters')
zlabel('depth of lake (meters)')
title(d);
d=d+1;
pause(0.005);
subplot(1,3,3);
plot(SA(:,2),SA(:,1),'.');
grid('on')
title('Exposed Shore Area')
xlabel('year')
ylabel('square kilometers (km^2)')
axis([0 50 0 200])
%%%%%%%%%%%%%%%%%%%%%%%%%% Lake Calculations %%%%%%%%%%%%%%%%%%%%%%%%%%%%
out=0;
LakeArea=0;
while m<a % Calculates the Lake water Area by number of vertices below zero
if(LakeCord(m,3)< -0.0000000001)
LakeArea=LakeArea+1;
end
m = m + 1;
end
m=1;
LakeArea=LakeArea*(1000000/6);
ground_out=ground_outO*(LakeArea/LakeAreaO);
out=evap*LakeArea+ground_out;
Salt=Salt+Salt_in;
volume = -(1000000/6)*sum(LakeCord(:,3));
while p < a % Calculates drop in water height of each vertical element
if (LakeCord(p,3) <= -0.0000000001)
LakeCord(p,3) = LakeCord(p,3) - ((input-out)/LakeArea)-Vcor/LakeArea;
end
Vcor=0;
p= p + 1;
end
p=1;
q=1;
while q < a % Corrects the volume if the change in volume falls below an element.
if LakeCord(q,3) > 0
Vcor = Vcor + (1000000/6)*LakeCord(q,3);
end
q=q+1;
end
subLakeCord=LakeCord(:,3); % Assigns 0 to all non negative depth values to the lake to calculate new area
subLakeCord(subLakeCord>0)=0;
LakeCord(:,3)=subLakeCord;
LakeCord;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Lake_Depth= -min(LakeCord(:,3));
Percent_of_Original_Volume = (volume/volumeO)*100;
acer_feet=0.00081071*volume;
LakeArea;
y;
i=i+1;
input;
out;
Salinity=Salt/volume;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Program Outputs%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Oringinal_Area_of_Lake_square_kilometers=(LakeAreaO/1000000)
exposed_Shoreline_Square_kilometers =(LakeAreaO-LakeArea)/1000000% Original Area of lake
Original_Acer_feet_of_Lake=volumeO/1233.48184 %Original Volume of lake (2006)
Final_Acer_Feet=acer_feet
volume; % volume of Lake
Final_Lake_Area_square_kilometers=LakeArea/1000000 % Exposed surface area of lake
Original_Lake_Depth_meters=Lake_DepthO
Lake_Depth_meters=Lake_Depth % Lake Depth and the end of the Simulation
Percent_of_Original_Volume
acer_feet;
Page updated
Google Sites
Report abuse