%% Matlab Livescript for processing DAS strainrate and water level collected from the NGLd well at Curtin University on April 8, 2022, Test 28.  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("HeadTableApr8.mat"); % Water level table data from logger in well
load('slug28_1Hz.mat'); %DAS matrix data decimated to 1 Hz (dec1Hz) for Test 28
HeadTable.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

% Compute displacement by integrating displacement rate, using original sampling rate for dt to get correct magnitude

DispRate=dec1Hz*fs;
Disp=double(cumtrapz(DispRate))*dt;
StrainRate=DispRate/10;                 %calculate strain rate by dividing by 10 m gauge length
Strain=Disp/10;                         %calculate strain by dividing by 10 m gauge length
MStrainRate=mean(DispRate(:,600:899),2);                %compute mean strain rate over channels
cmStrainRate=StrainRate-mean(StrainRate(:,50:400),2);   %remove common mode from Strain Rate
MStrain=mean(Strain(:,600:899),2);                      %compute mean strain over channels
cmStrain=Strain-mean(Strain(:,300:500),2);              %remove common mode from Strain

% Coordinate times between DAS data and water level change in the well
% Test 28 (Periodic slug) started on April 8, 2022 04:44 UTC

DASFStart28= datetime(2022,04,08,04,39,30,00,'TimeZone','UTC');   %time of DAS file start Test 28
TestStart28=datetime(2022,04,08,04,44,00,00,'TimeZone','UTC');    %time of well level data start Test 28
Tdas=DASFStart28+seconds(i*dt);
Tdas=Tdas';
Tdas.TimeZone=          'Asia/Shanghai';

%
%% Create Figure 2:
% Figure 2  Comparison of strain in more responsive (red) and less responsive (blue) strata. Strain rate is significantly different in each layer as the layers are expected to have marked differences in poro-elastic properties.

%Define colors for consistent presentation in red and blue
rbcolors = [1 0 0 
            0 0 1
            0 1 0 ];

%Define display time limit
xloc1=DASFStart28+minutes(4.5);
xloc2=DASFStart28+minutes(9);

%pick channels to compare
UpChan=694;  %712, %855  %662
LowChan=827;%790; %751, %867  %693
legup= num2str(round(D(UpChan))+" m");
leglow=num2str(round(D(LowChan))+" m");

f7=figure(2);  %plot strain in two sectons  This is Fig 2 in GRL
subplot(2,1,1)
    colororder(rbcolors) 
    v = pcolor(Tdas, D, cmStrain');
    set(v, 'EdgeColor', 'none')
    set(gca, 'clim', [-500 500]);
    colormap('jet');
    c7=colorbar; c7.Location="north";
    c7.Ruler.TickLabelFormat='%g nm/m'; c7.Limits=[-500,500];
    grid on; set(gca,'layer','top');
    %title("Strain Test 28")
    ylabel('Depth (m)')
    xlim([xloc1, xloc2])
    ylim([550,892])
    axis ij
    hold on;
        line([xloc1, xloc2], [D(UpChan),D(UpChan)], 'Color', 'r','LineWidth',2);
        line([xloc1, xloc2], [D(LowChan),D(LowChan)], 'Color', 'b','LineWidth',2);
        text(xloc1+seconds(1), D(UpChan)-25, legup  ) ;
        text(xloc1+seconds(1), D(LowChan)-25,leglow );
    hold off
subplot(2,1,2)
    colororder(rbcolors) 
    yyaxis left    
        plot(Tdas,detrend(cmStrain(:,UpChan)),'LineWidth',2) ; hold on
        ylabel('Strain (nanostrain)')
        ylim([-500,1000])
    yyaxis right
        plot(Tdas,detrend(cmStrain(:,LowChan)),'LineWidth',2); hold off
        legend(legup,leglow)
        ylabel('Strain (nanostrain)')
        ylim([-50,100])

    grid on; grid minor;
    xlim([DASFStart28+minutes(4.5), DASFStart28+minutes(9)])

%
%% Create Figure 3:
% Figure 3.  Strain rate in responsive zone that are correlated (red) and anticorrelated (blue) with change in water level measured in the well (black curves referenced to the right-hand water level axis) This test shows excellent repeatability through cycles and significant differences in the DAS strain rate recovered from the two depths.

%pick channels to compare
UpChan=712;  %712, %855  %662
LowChan=751; %751, %867  %693
legup= num2str(round(D(UpChan))+" m");
leglow=num2str(round(D(LowChan))+" m");
    
f3=figure(3);  %plot strain rate this is Figure 3 in GRL paper
subplot51=subplot(2,1,1);
    v = pcolor(Tdas, D, cmStrainRate');
    set(v, 'EdgeColor', 'none')
    set(gca, 'clim', [-50 50]);
    colormap('jet');
    c3=colorbar; c3.Location="north";
    c3.Ruler.TickLabelFormat='%g nm/m/s'; c3.Limits=[-50,50];
    grid on; set(gca,'layer','top');
    ylabel('Depth (m)')
    xlim([xloc1, xloc2])
    ylim([600,800])
    axis ij
    hold on;
        line([xloc1, xloc2], [D(UpChan),D(UpChan)], 'Color', 'r','LineWidth',2);
        line([xloc1, xloc2], [D(LowChan),D(LowChan)], 'Color', 'b','LineWidth',2);
        text(xloc1+seconds(1), D(UpChan)-10, legup  ) ;
        text(xloc1+seconds(1), D(LowChan)-10,leglow );
    hold off
subplot52=subplot(2,1,2);
    yyaxis left    
    plot(Tdas,detrend(cmStrainRate(:,UpChan)),'-r','LineWidth',2) ; hold on
    plot(Tdas,detrend(cmStrainRate(:,LowChan)),'-b','LineWidth',2); hold off
    ylabel('Strain Rate (nanostrain/s)')
    ylim([-300,300])
    yyaxis right
    plot(HeadTableApr8.DT,HeadTableApr8.H_m,'-k','LineWidth',2); %hold on;
    legend(legup,leglow,'Water Level')
    ylabel('Water Level (m)');
    xlabel('Minutes');
    set(subplot52,'YColor',[0 0 0]);
    grid on; grid minor;
    ylim([-1.5,1.5])
    xlim([DASFStart28+minutes(4.5), DASFStart28+minutes(9)])
%pick channels to compare
UpperDepth=732;
LowerDepth=772;
UpChan=round(UpperDepth/1.03+1);  
LowChan=round(LowerDepth/1.03+1); 