StateTL - updates of StateTL code for B&C calibration work

[kelleythompson]

A number of revisions have been made to the StateTL model engine code and input files prior to being used in the actual calibration work being undertaken by Brown and Caldwell. This issue will document those numerous revisions that were addressing a variety of issues and enhancements.

[kelleythompson]

Getting back to this...

Again, still had some problem with j349 producing spikes with step type releases. The following shows spiking in a very short channel due to bank storage compared to case where bank storage was turned off.

This wouldn't be a problem in calibration as the steps don't occur in full river flows. However, this does become a problem when doing release modeling either at low flows (when using the initial method with j349 of subtracting releases from river flows) or using a new method (programmed in the matlab code) to just model the flows by themselves but with full river celerity/dispersion (seems to be the only way to fully eliminate spike problem). So as we probably needed to do this for CoW implementation, thought needed to have this in place for calibration. Anyhow, so the biggest part of this work was to implement the j349 bank storage routine as much as possible into the matlab code.

Despite trying hard to do so, the exact j349 bank-storage algorithm could not be directly replicated.
The goal was to replicate the bank-storage aquifer timing response of J349 as closely as possible. The j349 response curve could be directly replicated, but the scaling of that curve to a volume as a function of change in stage could not be replicated (units were wrong and method created spikes etc etc). Therefore, the first ordinate of the glover solution to aquifer response was used to scale the j349 response curve, and this did replicate bank storage responses quite closely. One additional feature is to also cut off the response tail after it reaches 99% which significantly reduces model run times. General code:

T=transmissivity; % think in ft2/day SY=storagecoefficient; %unitless XL=aquiferwidth; %feet AL=alluviumlength; %miles rhours=1 %1 hour timestep

%j349 response curve ALPHA=(T / 24.) rhours / SY; %ALPHA as written in J349 (in wrong units) for NT=1:rsteps TIME=NT-.5; N=0; D=0; X1=1; %%%%%%%%%%%%%%%% % J349 Case 1 % DUSRF(NT)=-1 / sqrt (pi() ALPHA * TIME);

%%%%%%%%%%%%%%%%
% J349 Case 2 – bounded aquifer case
while X1>0.001
    N=N+1;
    DD=D;
    C1=(2*N-1)*pi()/(2.*XL);
    D=D+exp(-C1^2 * ALPHA * TIME);
    X1=abs(DD-D);
    DUSRF(NT)=(2./XL)*D;
    if N>100
        break
    end
end

end

% Glover bank storage F=H L T / sqrt(pi T/SS t) – infinite If bounded aquifer involves (1-exp(-(2XL)2/(4 T/SS t))with additional image wells, but the first ordinate is the same. %OK - only way could figure out with correct units was to unitize first ordinate and use first ordinate of glover solution to scale

T=T / (rhours 3600 24); %transmissibility of bank materials in ft2/s gloveralpha=2 AL 5280 T ./ sqrt(pi() T/SY * 1800); %1800=t(sec) at 0.5hr

bankurf=DUSRF/DUSRF(1); %unitize j349 curve bankurf=bankurf*gloveralpha; %j349 curve scaled to glover alpha

% cutting off end of tail urfwrapback=.99; %if wrapping back, wrap back portion bankurfcum=cumsum(bankurf); %ft2/s cburf=bankurfcum/bankurfcum(end); urfwrapids=find(cburf>=urfwrapback); if isempty(urfwrapids) urfwrapid=length(bankurf); else urfwrapid=urfwrapids(1); %ie includes first one over limit end bankurf=bankurf(1:urfwrapid);

Figure - Example of difference between j349 response (blue) and both glover infinite (yellow) and glover bounded (brown) solutions. This is the reason that the j349 curve was used although it had to be scaled with the first glover ordinate.

Figure - Example of responses over spinup time for j349 (brown) and new method (purple)

Figure - Example comparing response of j349 and new method without (light lines) and with evaporation (heavier lines). It is unsure if this figure is for that final implementation of new method or a slightly older version