%% Matlab Livescript for processing DAS strainrate and water level collected from the NGLd well at Curtin University on April 6, 2022 (Test 1).  This script supports the article "Distributed Water Storage Measured by Fiber Optic Distributed Acoustic Sensing" submitted to Geophysical Research Letters by Becker, Harris,Pevzner
%% This script was prepared using Matlab version 2022a.
%
% Load Data and set time zone.  Note that the original DAS data were collected at a sampling rate of 1 kHz.  The individual files were concatenated into a single matrix then decimated using the decimate command in matlab by a factor of 1000 to produce the 1 Hz sampling rate.

clear all; close all;
load("HeadTableApr6.mat"); % Water level table data from logger in well
load('slug1_4_1Hz.mat'); %DAS matrix data decimated to 1 Hz (dec1Hz) for Test 1
HeadTableApr6.DT.TimeZone =   'Asia/Shanghai';

% Set DAS sampling parameters

fs=1000;    %original sampling rate in Hz
dfs=1;      %decimated sampling rate in Hz
dt=1/dfs;   %time increment in seconds

i=1:size(dec1Hz,1);     %Count through time
j=1:size(dec1Hz,2);     %Count through channels
D=(j-1)*1.03;           %convert DAS measurement channel to depth; spatial measurement resolution = 1.03 m. 

% Compute displacement by integrating displacement rate, using original sampling rate for dt

DispRate=dec1Hz*fs;
Disp=double(cumtrapz(DispRate))*dt;

% Integrate displacement rate from 653 to 899 to capture total displacement over screened interval, then convert to displacement.

TotDispRate=sum(DispRate,2);
TotDisp=detrend(double(cumtrapz(TotDispRate))*dt);

% Coordinate times between DAS data and water level change in the well


DASFStart1= datetime(2022,04,06,02,39,00,00,'TimeZone','UTC');   %time of DAS file start Test 1
TestStart1=datetime(2022,04,06,02,40,00,00,'TimeZone','UTC');    %time of well level data start Test 1
Tdas=DASFStart1+seconds(i*dt);
Tdas=Tdas';
Tdas.TimeZone=          'Asia/Shanghai';

% Pick two channels/depths to compare in plots

%pick channels to compare
UpperDepth=732;
LowerDepth=772;
UpChan=round(UpperDepth/1.03+1);  
LowChan=round(LowerDepth/1.03+1); 

%

figure(1)
subplot(10,1,1:7)
    plot(Tdas,TotDisp,'LineWidth',2) ; 
    title('Total Displacement over Screened Interval')
    grid on; set(gca,'layer','top');
    ylabel('Total Displacement (nm)')
    xlim([DASFStart1+seconds(30), DASFStart1+minutes(2.5)])
    %ylim([0, 892])
subplot(10,1,9:10)
    plot(HeadTableApr6.DT,HeadTableApr6.H_m,'-k','LineWidth',2); 
    title('Measured Water Level in Well')
    ylabel('Water Level (m)')
    xlim([DASFStart1+seconds(30), DASFStart1+minutes(2.5)])

% Create plot showing displacmenet over the entire screened interval relative to head in the well.

xlim1=floor(4*60);xlim2=floor(12*60);
ylim1=floor(LowChan*.95);ylim2=floor(UpChan*1.05); %set limits according to channel
tplot=xlim1:xlim2; zplot=D(ylim1:ylim2);  % calc time(min) depth(m) for plots

%locate position in array at start time
WellIdx=find(HeadTableApr6.DT>TestStart1,120);   % locate first 510 seconds after test start
DASIdx=find(Tdas>DASFStart1,120);            % locate first 510 seconds after test start

% Strain and strain rate is calculated from displacement and displacement rate by dividing by the gauge length. In these experiments the gauge length was 10 m.

StrainRate=DispRate/10;
Strain=Disp/10;

% Create a waterfall plot of displacement rate to view overall DAS response.

%pick channels to compare
UpChan=855;  %712, %855
LowChan=867; %751, %867
    
xlim1=1;xlim2=floor(12*60);ylim1=floor(LowChan*.95);ylim2=floor(UpChan*1.05); %set limits according to channel
tplot=xlim1:xlim2; zplot=D(ylim1:ylim2);  % calc time(min) depth(m) for plots

figure(2)
subplot(10,1,1:7)
    v = pcolor(Tdas, D, StrainRate');
    set(v, 'EdgeColor', 'none')
    set(gca, 'clim', [-50 50]);
    colormap('jet');
    c7=colorbar; c7.Location="north";
    c7.Ruler.TickLabelFormat='%g nm/m/s'; c7.Limits=[-50,50];
    grid on; set(gca,'layer','top');
    ylabel('Depth (m)')
    xlim([DASFStart1+seconds(30), DASFStart1+minutes(2.5)])
    ylim([0, 892])
    axis ij
    hold on;
        line([DASFStart1+seconds(30), DASFStart1+minutes(3)], [653,653], 'Color', 'k','LineWidth',2,'LineStyle','--');
        text(DASFStart1+seconds(32),630,'Top of Screen')
    hold off
subplot(10,1,9:10)
    plot(HeadTableApr6.DT,HeadTableApr6.H_m,'-k','LineWidth',2); 
    ylabel('Water Level (m)')
    xlim([DASFStart1+seconds(30), DASFStart1+minutes(2.5)])

%