templatecem.m

Stanford的SRB实验室Quantitative Seismic Interpretation的免费MATLAB程序

%TEMPLATECEM RPT (rock physics template) for cemented sand
%        defans=templatecem(defans)
%
% Creates an cemented sand Rock Physics template plot as Vp/Vs vs. AI.
% This version is a combination of the soft sediment model for shale
% and the constant cement model for the sand.
% Output parameter defans records chosen parameters which can be used as
% defaults for next call.  Easiest way to start is call with no inputs and answer the
% prompts.
% Input DEFANS is optional.

%Written by Gary Mavko 02/04

% Mineral component elastic moduli
KQuartz=36.6e9; GQuartz=45.0e9; KClay=25e9; GClay=9e9; KFeldspar=75.6e9; GFeldspar=25.6e9;
KCalcite=76.8e9; GCalcite=32.0e9;
% Mineral component densities
RhoQuartz=2650; RhoClay=2580; RhoFeldspar=2630; RhoCalcite=2710;

    prompt = {'Effective Pressure (MPa)', ...
             'Fraction cement', ...
             'Quartz Fraction in sand; in shale', ...
             'Clay Fraction in sand; in shale', ...
             'Calcite Fraction in sand; in shale', ...
             'Feldspar Fraction in sand; in shale', ...
             'Shear Reduction Factor', ...
             'Brine Bulk Modulus',...
             'Brine Density', ...
             'Hydrocarbon Bulk Modulus',...
             'Hydrocarbon Density', ...
             'Clay Bulk Modulus',...
             'Clay Shear Modulus', ...
             'Clay Bulk density'};    
    if nargin<1,
        defans = {'40', ...
             '.03',...
             ['1';'0'],...
             ['0';'1'],...
             ['0';'0'],...
             ['0';'0'], ...
             '.5', ...
             '2.5e9',...
             '1020', ...
             '0.2e9',...
             '200', ...
             num2str(KClay), ...
             num2str(GClay), ...
             num2str(RhoClay)};         
     end;
 
    answer    = inputdlg(prompt,'Constant Cement Model',[1 1 2 2 2 2 1 1 1 1 1 1 1 1],defans);
    Pressure  = str2num(answer{1});
    fcement   = str2num(answer{2});
    xQuartz   = str2num(answer{3}); sQuartz   = xQuartz(1);   cQuartz   = xQuartz(2);
    xClay     = str2num(answer{4}); sClay     = xClay(1);     cClay     = xClay(2);
    xCalcite  = str2num(answer{5}); sCalcite  = xCalcite(1);  cCalcite  = xCalcite(2);
    xFeldspar = str2num(answer{6}); sFeldspar = xFeldspar(1); cFeldspar = xFeldspar(2);
    Sfact     = str2num(answer{7});
    Kbrine    = str2num(answer{8});
    Rhobrine  = str2num(answer{9});
    Kgas      = str2num(answer{10});
    Rhogas    = str2num(answer{11});
    KClay     = str2num(answer{12});
    GClay     = str2num(answer{13});
    RhoClay   = str2num(answer{14});
defans = answer;

Coord = 6;
colr = [0.6 0 0];

% First, compute soft shale properties
PhiC = .70;    % assumed critical porosity
phiplot = [.10 .20 .30 .40 .50 .60 .70];
[Phi,VpDry,VsDry,RHOBDry,VpSat,VsSat,RHOBSat]=softsediments(cClay,cFeldspar,cCalcite,Pressure,PhiC,Coord,Kbrine,Rhobrine,KClay,GClay,RhoClay,Sfact);
vpplot = interp1(Phi,VpSat,phiplot);
vsplot = interp1(Phi,VsSat,phiplot);
rhoplot = interp1(Phi,RHOBSat,phiplot);
h=plot(rhoplot.*vpplot,vpplot./vsplot,'o-'); set(h,'markersize',10);set(h,'color',colr); set(h,'markeredgecolor',colr);set(h,'markerfacecolor','w');
dx = .03*(max(vpplot.*rhoplot)-min(vpplot.*rhoplot));
dy = .01*(max(vpplot./vsplot)-min(vpplot./vsplot));
for k=1:length(phiplot)-1,h=text(vpplot(k)*rhoplot(k)+dx,vpplot(k)/vsplot(k)+dy,num2str(100*phiplot(k)));set(h,'color',colr); end;
h=text(vpplot(end)*rhoplot(end)+dx,vpplot(end)/vsplot(end)+dy,['Shale Por: ' num2str(100*phiplot(end))]);set(h,'color',colr);

hold on;

% Compute brine sand properties
PhiC = .40;    % assumed critical porosity
phiplot = [.10 .20 .30 .35];
[x,x,x,x,x,x,x,Phi,x,x,x,RHOBSat,VpSat,VsSat]=cementedsand(sClay,sFeldspar,sCalcite,PhiC,fcement,Kbrine,Rhobrine,Sfact);
vpplot = interp1(Phi,VpSat,phiplot);
vsplot = interp1(Phi,VsSat,phiplot);
rhoplot = interp1(Phi,RHOBSat,phiplot);
h=plot(rhoplot.*vpplot,vpplot./vsplot,'o-'); set(h,'markersize',10);set(h,'color',colr);set(h,'markeredgecolor',colr);set(h,'markerfacecolor','w');
for k=1:length(phiplot)-1,h=text(vpplot(k)*rhoplot(k)+dx,vpplot(k)/vsplot(k)+dy,num2str(100*phiplot(k)));set(h,'color',colr);end;
h=text(vpplot(end)*rhoplot(end)+dx,vpplot(end)/vsplot(end)+dy,['Sand Por: ' num2str(100*phiplot(end))]);set(h,'color',colr);

% Repeat sand model, looping over gas saturations
Sw = [1 .9 .8 .7 .6 .5 .4 .3 .2 .1 0];
PhiC = .40;    % assumed critical porosity
phiplot = [.10 .20 .30 .35];
for j = 1:length(Sw),
    Kfl   = 1/(Sw(j)/Kbrine + (1-Sw(j))/Kgas);
    Rhofl = Sw(j)*Rhobrine + (1-Sw(j))*Rhogas;
    [x,x,x,x,x,x,x,Phi,x,x,x,RHOBSat,VpSat,VsSat]=cementedsand(sClay,sFeldspar,sCalcite,PhiC,fcement,Kfl,Rhofl,Sfact);
    vpplot(j,:)  = interp1(Phi,VpSat,phiplot);
    vsplot(j,:)  = interp1(Phi,VsSat,phiplot);
    rhoplot(j,:) = interp1(Phi,RHOBSat,phiplot);
end;

for k = 1:length(phiplot),
    h=plot(rhoplot(:,k).*vpplot(:,k),vpplot(:,k)./vsplot(:,k),'o-'); set(h,'markersize',10);set(h,'color',colr);set(h,'markeredgecolor',colr);set(h,'markerfacecolor','w');
end;

h=text(rhoplot(1,end).*vpplot(1,end) -dx ,vpplot(1,end)./vsplot(1,end),'Sw= 1');set(h,'horizontalalignment','right');set(h,'color',colr);
h=text(rhoplot(end,end).*vpplot(end,end) -dx ,vpplot(end,end)./vsplot(end,end),'Sw= 0');set(h,'horizontalalignment','right');set(h,'color',colr);

Ks=0.5.*(sQuartz.*KQuartz+sFeldspar.*KFeldspar+sCalcite.*KCalcite+sClay.*KClay+1./(sQuartz./KQuartz+sFeldspar./KFeldspar+sCalcite./KCalcite+sClay./KClay));
Gs=0.5.*(sQuartz.*GQuartz+sFeldspar.*GFeldspar+sCalcite.*GCalcite+sClay.*GClay+1./(sQuartz./GQuartz+sFeldspar./GFeldspar+sCalcite./GCalcite+sClay./GClay));
RHOs=sQuartz.*RhoQuartz+sFeldspar.*RhoFeldspar+sCalcite.*RhoCalcite+sClay.*RhoClay;

Kc=0.5.*(cQuartz.*KQuartz+cFeldspar.*KFeldspar+cCalcite.*KCalcite+cClay.*KClay+1./(cQuartz./KQuartz+cFeldspar./KFeldspar+cCalcite./KCalcite+cClay./KClay));
Gc=0.5.*(cQuartz.*GQuartz+cFeldspar.*GFeldspar+cCalcite.*GCalcite+cClay.*GClay+1./(cQuartz./GQuartz+cFeldspar./GFeldspar+cCalcite./GCalcite+cClay./GClay));
RHOc=cQuartz.*RhoQuartz+cFeldspar.*RhoFeldspar+cCalcite.*RhoCalcite+cClay.*RhoClay;

f = 1e-6;
g = 1e-9;
params = {['Cement Model']; ...
         ['Sand: Qz=',num2str(sQuartz),'; C=',num2str(sClay),'; F=',num2str(sFeldspar),'; Ca=',num2str(sCalcite)]; ...
         ['Shale: Qz=',num2str(cQuartz),': C=',num2str(cClay),'; F=',num2str(cFeldspar),'; Ca=',num2str(cCalcite)]; ...
         ['Ksand=',num2str(g*Ks),' GPa']; ...
         ['Gsand=',num2str(g*Gs),' GPa']; ...
         ['Rhosand=',num2str(RHOs),' Kg/m3']; ...
         ['']; ...
         ['Kshale=',num2str(g*Kc),' GPa']; ...
         ['Gshale=',num2str(g*Gc),' GPa']; ...
         ['Rhoshale=',num2str(RHOc),' Kg/m3']; ...
         ['']; ...
         ['Khydrocarb=',num2str(g*Kgas),' GPa']; ...
         ['Rhohydrocarb=',num2str(Rhogas),' Kg/m3']; ...
         ['']; ...
         ['Kbrine=',num2str(g*Kbrine),' GPa']; ...
         ['Rhobrine=',num2str(Rhobrine),' Kg/m3']; ...
         ['']; ...
         ['frac cement=',num2str(fcement)]; ...
         ['Press=',num2str(Pressure),' MPa']};
         
 a = axis; dxx = 0.1*(a(2)-a(1));dyy = 0.1*(a(4)-a(3));
 axis([a(1)-dxx,a(2)+dxx,a(3)-dyy,a(4)+dyy]);
 a = axis; dxx = 0.1*(a(2)-a(1));dyy = 0.1*(a(4)-a(3));
 h=text(a(2)-2.7*dxx,a(4),params);set(h,'fontsize',8,'horizontalalignment','left','verticalalignment','top');set(h,'color',colr);
 
 xlabel('Acoustic Impedance');
 ylabel('Vp/Vs');