diff --git a/.gitignore b/.gitignore
index db01465..6caddde 100644
--- a/.gitignore
+++ b/.gitignore
@@ -18,3 +18,4 @@ Projects/1_LabData/1_Attenuation/not_used/
*.m~
*.idea
.vbr_test_data_dir
+*.zip
\ No newline at end of file
diff --git a/Projects/1_LabData/1_Attenuation/FitData_YT16.m b/Projects/1_LabData/1_Attenuation/FitData_YT16.m
index ea112bb..dcfb73b 100644
--- a/Projects/1_LabData/1_Attenuation/FitData_YT16.m
+++ b/Projects/1_LabData/1_Attenuation/FitData_YT16.m
@@ -18,6 +18,7 @@ function FitData_YT16()
path_to_top_level_vbr='../../../';
addpath(path_to_top_level_vbr)
vbr_init
+ addpath('./functions')
plot_visc();
plot_Q();
@@ -32,7 +33,6 @@ function plot_Q()
viscData = loadYT2016visc();
Qdata=loadYT2016Q();
- %figure('DefaultAxesFontSize',12,'PaperPosition',[0,0,6,2.5],'PaperPositionMode','manual')
fig=figure('Position', [10 10 600 300],'PaperPosition',[0,0,6,2.5],'PaperPositionMode','manual');
OutVBR=struct();
if viscData.has_data && Qdata.Qinv.has_data
@@ -63,7 +63,7 @@ function plot_Q()
[T_Cvisc,I]=sort(T_Cvisc); eta=eta(I);
% set this sample's viscosity parameters
- VBR.in.viscous.xfit_premelt=setBorneolParams();
+ VBR.in.viscous.xfit_premelt=setBorneolViscParams();
VBR.in.viscous.xfit_premelt.dg_um_r=dg;
VBR.in.viscous.xfit_premelt.Tr_K=T_Cvisc(1)+273;
VBR.in.viscous.xfit_premelt.eta_r=eta(1);
@@ -300,56 +300,3 @@ function plot_visc()
end
end
-function params = setBorneolParams()
- % set the general viscous parameters for borneol.
- % near-solidus and melt effects
- params.alpha=25;
- params.T_eta=0.94; % eqn 17,18- T at which homologous T for premelting.
- params.gamma=5;
- % flow law constants for YT2016
- params.Tr_K=23+273; % reference temp [K]
- params.Pr_Pa=0; % reference pressure [Pa]
- params.eta_r=7e13;% reference eta (eta at Tr_K, Pr_Pa)
- params.H=147*1e3; % activation energy [J/mol]
- params.V=0; % activation vol [m3/mol]
- params.R=8.314; % gas constant [J/mol/K]
- params.m=2.56; % grain size exponent
- params.dg_um_r=34.2 ; % eference grain size [um]
-end
-
-function [G,dGdT,dGdT_ave] = YT16_E(T_C)
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % [G,dGdT,dGdT_ave]= YT16_E(T_C)
- %
- % YT2016 Figure 6 caption polynomial fit of E vs T, T in degrees C
- %
- % parameters
- % ----------
- % T_C temperature in deg C, any size array
- %
- % output
- % ------
- % G modulus, GPa. same size as T_C
- % dGdT temperature derivative of G at T_C, same size as T_C. [GPa/C]
- % dGdT_ave average temp derivative over range of T_C given, scalar. [GPa/C].
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
- a0=2.5943;
- a1=2.6911*1e-3;
- a2=-2.9636*1e-4;
- a3 =1.4932*1e-5;
- a4 =-2.9351*1e-7;
- a5 =1.8997*1e-9;
- a = [a5,a4,a3,a2,a1,a0];
-
- G=zeros(size(T_C));
- dGdT=zeros(size(T_C));
- for iTC=1:numel(T_C)
- Gpoly = polyval(a,T_C(iTC));
- dGdTpoly = polyval(a(1:end-1),T_C(iTC));
- G(iTC)=sum(Gpoly(:));
- dGdT(iTC)=sum(dGdTpoly(:));
- end
- dGdT_ave=mean(dGdT);
-
-end
diff --git a/Projects/1_LabData/1_Attenuation/FitData_YT24.m b/Projects/1_LabData/1_Attenuation/FitData_YT24.m
new file mode 100644
index 0000000..d0a9461
--- /dev/null
+++ b/Projects/1_LabData/1_Attenuation/FitData_YT24.m
@@ -0,0 +1,214 @@
+function VBRc_results = FitData_YT24()
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % FitData_YT24()
+ %
+ % Plots viscosity, modulus, Qinv for borneol near and above the solidus
+ % temperature following premelting scaling of Yamauchi and Takei, 2024, JGR.
+ %
+ % Parameters
+ % ----------
+ % use_data
+ % if 1, will attempt to fetch data from zenondo
+ %
+ % Output
+ % ------
+ % figures to screen and to Projects/1_LabData/1_Attenuation/figures/
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ % put VBR in the path
+ vbr_path = getenv('vbrdir');
+ if isempty(vbr_path)
+ vbr_path='../../../';
+ end
+ addpath(vbr_path)
+ vbr_init
+
+ addpath('./functions')
+
+ % check on data, download if needed
+ full_data_dir = download_YT24('data');
+
+ % get VBRc result and plot along the way
+ VBRc_results = plot_fig7(full_data_dir);
+
+end
+
+
+function VBRc_results = plot_fig7(full_data_dir);
+
+ figure('Position', [10 10 500 600],'PaperPosition',[0,0,7,7*6/5],'PaperPositionMode','manual');
+ yheight = 0.4;
+ xpos = 0.1;
+ ypos = 0.1;
+ xwidth = 0.6;
+ ax_Q = axes('position', [xpos, ypos, xwidth, yheight]);
+ ax_E = axes('position', [xpos, ypos + yheight, xwidth, yheight]);
+ leg_pos = [xpos+xwidth+0.05, ypos, 0.2, yheight];
+
+
+ load(fullfile('data','YT16','table3.mat'));
+ visc_data.table3_H=table3_H.table3_H;
+
+ combined_data = YT24_load_fig7_combined_data(full_data_dir);
+
+ n_exps = numel(combined_data);
+
+ freq_range = logspace(-6,9,100);
+ Tsol = 43.0;
+
+ % get VBRc results
+ VBRc_results = struct();
+ for i_exp = 1:n_exps
+ data = combined_data(i_exp);
+
+ T = data.T; % celcius
+ dg = data.dg_um;
+ phi = data.phi;
+ rho = data.rho_kgm3;
+
+ VBR.in.elastic.methods_list={'anharmonic';};
+ VBR.in.viscous.methods_list={'xfit_premelt'};
+ VBR.in.anelastic.methods_list={'xfit_premelt'};
+ VBR.in.anelastic.xfit_premelt.include_direct_melt_effect = 1.0;
+ % VBR.in.anelastic.xfit_premelt.tau_pp=2*1e-5;
+ % plotting with correction for poro-elastic effect applied, so
+ % set effect to 0 here.
+ VBR.in.anelastic.xfit_premelt.poro_Lambda = 0.0;
+
+ % set viscous params
+ VBR.in.viscous.xfit_premelt=setBorneolViscParams(); % YT2016 values, sample 41
+ VBR.in.viscous.xfit_premelt.alpha = 32; % YT2024 actually fit this
+
+ % set anharmonic conditions
+ VBR.in.elastic.anharmonic=Params_Elastic('anharmonic');
+ % extract reference modulus (E_normed = E / E_u)
+ [Gu_o,dGdT,dGdT_ave]= YT16_E(T);
+ Gu_o = Gu_o * 1e9;
+
+ % Gu_o is for a given T, specify at elevated TP to skip anharmonic scaling
+ VBR.in.elastic.Gu_TP = Gu_o;
+ VBR.in.elastic.quiet = 1; % not bothering with K, avoid printing the poisson warning
+
+ % set experimental conditions
+ VBR.in.SV.T_K = T+273 ;
+ sz=size(VBR.in.SV.T_K);
+ VBR.in.SV.dg_um= dg .* ones(sz); %dg.* ones(sz);
+ VBR.in.SV.P_GPa = 1.0132e-04 .* ones(sz); % pressure [GPa]
+ VBR.in.SV.rho =rho .* ones(sz); % density [kg m^-3]
+ VBR.in.SV.sig_MPa = 1000 .* ones(sz)./1e6; % differential stress [MPa]
+ VBR.in.SV.Tsolidus_K = Tsol + 273 ;
+
+ VBR.in.SV.phi = phi * ones(sz); % melt fraction
+ VBR.in.SV.Ch2o_0=zeros(sz);
+
+ VBR.in.SV.f=freq_range;
+ [VBR_bysamp] = VBR_spine(VBR);
+ results.Qinv=VBR_bysamp.out.anelastic.xfit_premelt.Qinv;
+ results.E=VBR_bysamp.out.anelastic.xfit_premelt.M;
+ results.f = VBR.in.SV.f;
+ results.f_normed = VBR_bysamp.out.anelastic.xfit_premelt.f_norm;
+ results.E_normed = results.E / Gu_o;
+ results.T = T;
+ results.Tn = VBR.in.SV.T_K / VBR.in.SV.Tsolidus_K ;
+ results.VBR = VBR_bysamp;
+ VBRc_results(i_exp) = results;
+ end
+
+
+ for i_exp = 1:n_exps
+ if i_exp == 2
+ linesty = '--';
+ else
+ linesty = '-';
+ end
+ results = VBRc_results(i_exp);
+ rgb = vbr_categorical_color(i_exp);
+ Tnlab = num2str(results.Tn, 3+(results.Tn>=1));
+ axes(ax_E)
+ hold_if(i_exp)
+ semilogx(results.f_normed, results.E_normed ,'color', rgb, ...
+ 'linewidth', 1.5, 'linestyle', linesty)
+ axes(ax_Q)
+ hold_if(i_exp)
+ loglog(results.f_normed, results.Qinv , 'color', rgb, ...
+ 'linewidth', 1.5, 'displayname', Tnlab, 'linestyle', linesty)
+ end
+
+ % plot experimental results
+ delta_systematic = 0.031;
+
+ for i_exp = 1:n_exps
+ data = combined_data(i_exp);
+ dporo = data.dporo * data.phi;
+ mod_fac = 1 / ( 1 - delta_systematic) / (1 - dporo); % correction factor
+ rgb = vbr_categorical_color(i_exp);
+ if data.Tn >= 1
+ mrkr = 'o';
+ mrkr_sz = 5;
+ mkr_fc = 'w';
+ else
+ mrkr = '.';
+ mrkr_sz=15;
+ mkr_fc = rgb;
+ end
+ Tnlab = num2str(data.Tn, 3+(data.Tn>=1));
+ axes(ax_E)
+ hold_if(i_exp)
+ semilogx(data.f_normed, data.E_normed * mod_fac, mrkr, ...
+ 'markersize', mrkr_sz, 'color', rgb, 'displayname', Tnlab, ...
+ 'MarkerFaceColor', mkr_fc)
+ axes(ax_Q)
+ hold_if(i_exp)
+ loglog(data.f_normed, data.Qinv, mrkr, 'markersize', mrkr_sz, ....
+ 'color', rgb, 'displayname', Tnlab, 'MarkerFaceColor', mkr_fc)
+ end
+ axes(ax_Q)
+ leg = legend('Location', 'eastoutside', 'Orientation', 'vertical', ...
+ 'AutoUpdate','off','title','T_n','NumColumns',1);
+
+ axes(ax_E)
+ xticklabels([])
+ xlim([1e-2,1e9])
+ ylim([0,1.1])
+ xlabel('f / f_M')
+ ylabel("E/E_u")
+
+ axes(ax_Q)
+ xlim([1e-2,1e9])
+ ylim([1e-3,2])
+ xlabel('f / f_M')
+ ylabel("Q^-^1")
+
+ set(leg, 'Position', leg_pos)
+
+ saveas(gcf,'./figures/YT24_MQ.png')
+end
+
+
+
+
+function unzipped_data_dir = download_YT24(datadir)
+ % downloads, unzips data from zenodo
+
+ unzipped_data_dir = fullfile(datadir,'YT24');
+ if exist(datadir) ~= 7
+ mkdir(datadir)
+ mkdir(unzipped_data_dir)
+ end
+
+ YT24data_url = 'https://zenodo.org/api/records/10201689/files-archive';
+ zipfilename = '10201689.zip';
+
+ local_data = fullfile(unzipped_data_dir, zipfilename);
+ if exist(local_data) ~= 2
+ urlwrite(YT24data_url, zipfilename)
+ movefile(zipfilename, unzipped_data_dir)
+ unzip(local_data, unzipped_data_dir)
+ end
+end
+
+function hold_if(iterval)
+ if iterval > 1
+ hold all
+ end
+end
\ No newline at end of file
diff --git a/Projects/1_LabData/1_Attenuation/data/YT16/table3.mat b/Projects/1_LabData/1_Attenuation/data/YT16/table3.mat
new file mode 100644
index 0000000..a6ce9c3
Binary files /dev/null and b/Projects/1_LabData/1_Attenuation/data/YT16/table3.mat differ
diff --git a/Projects/1_LabData/1_Attenuation/data/YT24/readme.pdf b/Projects/1_LabData/1_Attenuation/data/YT24/readme.pdf
new file mode 100644
index 0000000..5e2526d
Binary files /dev/null and b/Projects/1_LabData/1_Attenuation/data/YT24/readme.pdf differ
diff --git a/Projects/1_LabData/1_Attenuation/functions/YT16_E.m b/Projects/1_LabData/1_Attenuation/functions/YT16_E.m
new file mode 100644
index 0000000..fac0f05
--- /dev/null
+++ b/Projects/1_LabData/1_Attenuation/functions/YT16_E.m
@@ -0,0 +1,36 @@
+function [G,dGdT,dGdT_ave] = YT16_E(T_C)
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % [G,dGdT,dGdT_ave]= YT16_E(T_C)
+ %
+ % YT2016 Figure 6 caption polynomial fit of E vs T, T in degrees C
+ %
+ % parameters
+ % ----------
+ % T_C temperature in deg C, any size array
+ %
+ % output
+ % ------
+ % G modulus, GPa. same size as T_C
+ % dGdT temperature derivative of G at T_C, same size as T_C. [GPa/C]
+ % dGdT_ave average temp derivative over range of T_C given, scalar. [GPa/C].
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ a0=2.5943;
+ a1=2.6911*1e-3;
+ a2=-2.9636*1e-4;
+ a3 =1.4932*1e-5;
+ a4 =-2.9351*1e-7;
+ a5 =1.8997*1e-9;
+ a = [a5,a4,a3,a2,a1,a0];
+
+ G=zeros(size(T_C));
+ dGdT=zeros(size(T_C));
+ for iTC=1:numel(T_C)
+ Gpoly = polyval(a,T_C(iTC));
+ dGdTpoly = polyval(a(1:end-1),T_C(iTC));
+ G(iTC)=sum(Gpoly(:));
+ dGdT(iTC)=sum(dGdTpoly(:));
+ end
+ dGdT_ave=mean(dGdT);
+
+end
diff --git a/Projects/1_LabData/1_Attenuation/functions/YT24_load_fig7_combined_data.m b/Projects/1_LabData/1_Attenuation/functions/YT24_load_fig7_combined_data.m
new file mode 100644
index 0000000..43cc2f1
--- /dev/null
+++ b/Projects/1_LabData/1_Attenuation/functions/YT24_load_fig7_combined_data.m
@@ -0,0 +1,107 @@
+function combined_data_final = YT24_load_fig7_combined_data(full_data_dir)
+
+ combined_data = struct();
+ % pull all data for sample 41 first T cycle
+ sample_id= 41;
+ [current_run, sample_info] = load_sample_data(sample_id, full_data_dir);
+
+ start_end = sample_info.first_Tcycle_index_start_stop;
+ i_Tstart = start_end(1);
+ i_Tend = start_end(2);
+
+ i_data = 1;
+ Tsol = 43;
+ Tvals = [];
+ for iT = i_Tstart: i_Tend
+ combined_data(i_data) = extract_single_experiment(current_run, sample_info, iT, Tsol, sample_id);
+ Tvals = [Tvals, combined_data(i_data).Tn];
+ i_data = i_data + 1;
+ end
+
+ for sample_id = 42:8:50 % 42, 50
+ [current_run, sample_info] = load_sample_data(sample_id, full_data_dir);
+ start_end = sample_info.first_Tcycle_index_start_stop;
+ i_Tstart = start_end(1);
+ i_Tend = start_end(2);
+ for iT = i_Tstart: i_Tend
+ if current_run(iT).T >= Tsol
+ combined_data(i_data) = extract_single_experiment(current_run, sample_info, iT, Tsol, sample_id);
+ Tvals = [Tvals, combined_data(i_data).Tn];
+ i_data = i_data + 1;
+ end
+ end
+ end
+
+
+
+ % finally, re-order based on homologous temperature
+ [sorted_T, idx] = sort(Tvals);
+ for isamp = 1:numel(combined_data)
+ combined_data_final(isamp) = combined_data(idx(isamp));
+ end
+
+end
+
+function single_exp = extract_single_experiment(current_run, sample_info, iT, Tsol, sample_id)
+ T = current_run(iT).T; % celcius
+ cdata = current_run(iT).data;
+ single_exp.T = T;
+ single_exp.Tn = (T + 273.0) / (Tsol+273); % homologous temperature
+ single_exp.f = cdata(:,1);
+ single_exp.f_normed = cdata(:,2);
+ single_exp.E = cdata(:,3);
+ single_exp.E_normed = cdata(:,4);
+ single_exp.Qinv = cdata(:,5);
+
+
+ phi = sample_info.phi * (T >= Tsol);
+ single_exp.phi = phi;
+ single_exp.dg_um = sample_info.dg_um;
+ single_exp.sample_id = sample_id;
+ single_exp.rho_kgm3 = sample_info.rho_kgm3;
+ single_exp.dporo = (phi > 0) * sample_info.dporo;
+end
+
+function [run_data, sample_info] = load_sample_data(sample_number, full_data_dir)
+
+ sample_str = num2str(sample_number);
+ sample_file = fullfile(full_data_dir, ['anela',sample_str,'.mat']);
+ struct_name = ['anela',sample_str];
+ sample_data = load(sample_file); % e.g., sample_data.anela42
+ data = getfield(sample_data, struct_name); % e.g., anela42
+
+ run_data = data.run;
+ sample_info = get_sample_info(sample_number);
+end
+
+
+
+function sample_data = get_sample_info(sample_number)
+
+ if sample_number == 50
+ sample_data.dg_um = 47.9;
+ sample_data.rho_kgm3 = 1.0111;
+ sample_data.phi = 0.0368;
+ sample_data.eta_r = 1061*1e12;
+ sample_data.first_Tcycle_index_start_stop = [2, 9];
+ sample_data.nT = 8;
+ sample_data.dporo = 0.135;
+ elseif sample_number == 42
+ sample_data.dg_um = 46.3;
+ sample_data.rho_kgm3 = 1.0111;
+ sample_data.phi = 0.016;
+ sample_data.eta_r = 132 * 1e12;
+ sample_data.first_Tcycle_index_start_stop = [2, 8];
+ sample_data.nT = 7;
+ sample_data.dporo = 0.0585;
+ elseif sample_number == 41
+ sample_data.dg_um = 34.2;
+ sample_data.rho_kgm3 = 1.0111;
+ sample_data.phi = 0.004;
+ sample_data.eta_r = 1433 * 1e12;
+ sample_data.first_Tcycle_index_start_stop = [3, 10];
+ sample_data.nT = 8;
+ sample_data.dporo = 0.0178;
+ end
+
+end
\ No newline at end of file
diff --git a/Projects/1_LabData/1_Attenuation/functions/setBorneolViscParams.m b/Projects/1_LabData/1_Attenuation/functions/setBorneolViscParams.m
new file mode 100644
index 0000000..a1f39da
--- /dev/null
+++ b/Projects/1_LabData/1_Attenuation/functions/setBorneolViscParams.m
@@ -0,0 +1,16 @@
+function params = setBorneolViscParams()
+ % set the general viscous parameters for borneol for YT2016
+ % near-solidus and melt effects
+ params.alpha=25;
+ params.T_eta=0.94; % eqn 17,18- T at which homologous T for premelting.
+ params.gamma=5;
+ % flow law constants for YT2016
+ params.Tr_K=23+273; % reference temp [K]
+ params.Pr_Pa=0; % reference pressure [Pa]
+ params.eta_r=7e13;% reference eta (eta at Tr_K, Pr_Pa)
+ params.H=147*1e3; % activation energy [J/mol]
+ params.V=0; % activation vol [m3/mol]
+ params.R=8.314; % gas constant [J/mol/K]
+ params.m=2.56; % grain size exponent
+ params.dg_um_r=34.2 ; % eference grain size [um]
+end
\ No newline at end of file
diff --git a/Projects/vbr_core_examples/CB_014_xfit_premelt_extended.m b/Projects/vbr_core_examples/CB_014_xfit_premelt_extended.m
new file mode 100644
index 0000000..a9451ba
--- /dev/null
+++ b/Projects/vbr_core_examples/CB_014_xfit_premelt_extended.m
@@ -0,0 +1,74 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% CB_014_xfit_premelt_extended.m
+%
+% Calls VBR using xfit_premelt method including the direct melt effect,
+% from:
+% Hatsuki Yamauchi and Yasuko Takei, JGR 2024, "Effect of Melt on
+% Polycrystal Anelasticity", https://doi.org/10.1029/2023JB027738
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%% put VBR in the path %%
+ clear; close all
+ path_to_top_level_vbr='../../';
+ addpath(path_to_top_level_vbr)
+ vbr_init
+
+%% write method list %%
+ VBR.in.elastic.methods_list={'anharmonic'};
+ VBR.in.anelastic.methods_list={'xfit_premelt'};
+ VBR.in.anelastic.xfit_premelt.include_direct_melt_effect = 1;
+
+ % load anharmonic parameters, adjust Gu_0_ol and derivatives to match YT2016
+ VBR.in.elastic.anharmonic.Gu_0_ol=72.45; %[GPa]
+ VBR.in.elastic.anharmonic.dG_dT = -10.94*1e6; % Pa/C (equivalent ot Pa/K)
+ VBR.in.elastic.anharmonic.dG_dP = 1.987; % GPa / GPa
+
+
+ %% Define the Thermodynamic State %%
+ VBR.in.SV.T_K=1200:5:1500;
+ VBR.in.SV.T_K=VBR.in.SV.T_K+273;
+ sz=size(VBR.in.SV.T_K); % temperature [K]
+ VBR.in.SV.P_GPa = full_nd(2.5, sz); % pressure [GPa]
+
+ Tn_cases = [.96, .98, 1.0, 1.1, 1.1, 1.1, 1.1];
+ phi_cases = [0., 0., 0., 0.1, 1, 2, 4]*0.01;
+
+ % remaining state variables (ISV)
+ VBR.in.SV.dg_um=full_nd(.004*1e6, sz); % grain size [um]
+ VBR.in.SV.rho = full_nd(3300, sz); % density [kg m^-3]
+ VBR.in.SV.sig_MPa = full_nd(1, sz); % differential stress [MPa]
+ VBR.in.SV.f = 1; % 1 Hz
+
+
+ f1=figure('PaperPosition',[0,0,6,4],'PaperPositionMode','manual');
+ nTn = numel(Tn_cases);
+ VBR_results = struct();
+
+ for iTn = 1:nTn
+ VBRi = VBR;
+ VBRi.in.SV.phi = full_nd(phi_cases(iTn), sz);
+ VBRi.in.SV.Tsolidus_K = VBR.in.SV.T_K / Tn_cases(iTn);
+ [VBRi] = VBR_spine(VBRi) ;
+
+ results.Q = VBRi.out.anelastic.xfit_premelt.Q;
+ VBR_results(iTn) = results;
+
+ dname = [num2str(Tn_cases(iTn)), ', ', num2str(phi_cases(iTn))];
+ figure(f1)
+ if iTn > 1
+ hold all
+ end
+ plot(VBR.in.SV.T_K - 273, results.Q, 'displayname', dname, 'linewidth', 1.5)
+ legend('Location','eastoutside','title', '(Tn, phi)')
+ xlabel('Temperature [C]', 'fontsize', 12)
+ ylabel('Qs', 'fontsize', 12)
+ ylim([0, 200])
+
+ yticks(0:20:200)
+
+ end
+
+
+ saveas(gcf,'./figures/CB_014_xfit_premelt_extended.png')
+
+
diff --git a/docs/_pages/examples/CB_014_xfit_premelt_extended.md b/docs/_pages/examples/CB_014_xfit_premelt_extended.md
new file mode 100644
index 0000000..bf84e72
--- /dev/null
+++ b/docs/_pages/examples/CB_014_xfit_premelt_extended.md
@@ -0,0 +1,86 @@
+---
+permalink: /examples/CB_014_xfit_premelt_extended/
+title: ""
+---
+
+# CB_014_xfit_premelt_extended.m
+## output figures
+
+{:class="img-responsive"}
+## contents
+```matlab
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% CB_014_xfit_premelt_extended.m
+%
+% Calls VBR using xfit_premelt method including the direct melt effect,
+% from:
+% Hatsuki Yamauchi and Yasuko Takei, JGR 2024, "Effect of Melt on
+% Polycrystal Anelasticity", https://doi.org/10.1029/2023JB027738
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%% put VBR in the path %%
+ clear; close all
+ path_to_top_level_vbr='../../';
+ addpath(path_to_top_level_vbr)
+ vbr_init
+
+%% write method list %%
+ VBR.in.elastic.methods_list={'anharmonic'};
+ VBR.in.anelastic.methods_list={'xfit_premelt'};
+ VBR.in.anelastic.xfit_premelt.include_direct_melt_effect = 1;
+
+ % load anharmonic parameters, adjust Gu_0_ol and derivatives to match YT2016
+ VBR.in.elastic.anharmonic.Gu_0_ol=72.45; %[GPa]
+ VBR.in.elastic.anharmonic.dG_dT = -10.94*1e6; % Pa/C (equivalent ot Pa/K)
+ VBR.in.elastic.anharmonic.dG_dP = 1.987; % GPa / GPa
+
+
+ %% Define the Thermodynamic State %%
+ VBR.in.SV.T_K=1200:5:1500;
+ VBR.in.SV.T_K=VBR.in.SV.T_K+273;
+ sz=size(VBR.in.SV.T_K); % temperature [K]
+ VBR.in.SV.P_GPa = full_nd(2.5, sz); % pressure [GPa]
+
+ Tn_cases = [.96, .98, 1.0, 1.1, 1.1, 1.1, 1.1];
+ phi_cases = [0., 0., 0., 0.1, 1, 2, 4]*0.01;
+
+ % remaining state variables (ISV)
+ VBR.in.SV.dg_um=full_nd(.004*1e6, sz); % grain size [um]
+ VBR.in.SV.rho = full_nd(3300, sz); % density [kg m^-3]
+ VBR.in.SV.sig_MPa = full_nd(1, sz); % differential stress [MPa]
+ VBR.in.SV.f = 1; % 1 Hz
+
+
+ f1=figure('PaperPosition',[0,0,6,4],'PaperPositionMode','manual');
+ nTn = numel(Tn_cases);
+ VBR_results = struct();
+
+ for iTn = 1:nTn
+ VBRi = VBR;
+ VBRi.in.SV.phi = full_nd(phi_cases(iTn), sz);
+ VBRi.in.SV.Tsolidus_K = VBR.in.SV.T_K / Tn_cases(iTn);
+ [VBRi] = VBR_spine(VBRi) ;
+
+ results.Q = VBRi.out.anelastic.xfit_premelt.Q;
+ VBR_results(iTn) = results;
+
+ dname = [num2str(Tn_cases(iTn)), ', ', num2str(phi_cases(iTn))];
+ figure(f1)
+ if iTn > 1
+ hold all
+ end
+ plot(VBR.in.SV.T_K - 273, results.Q, 'displayname', dname, 'linewidth', 1.5)
+ legend('Location','eastoutside','title', '(Tn, phi)')
+ xlabel('Temperature [C]', 'fontsize', 12)
+ ylabel('Qs', 'fontsize', 12)
+ ylim([0, 200])
+
+ yticks(0:20:200)
+
+ end
+
+
+ saveas(gcf,'./figures/CB_014_xfit_premelt_extended.png')
+
+
+```
diff --git a/docs/_pages/examples/vbrcore.md b/docs/_pages/examples/vbrcore.md
index d2e9e27..03b44e0 100644
--- a/docs/_pages/examples/vbrcore.md
+++ b/docs/_pages/examples/vbrcore.md
@@ -18,3 +18,4 @@ Simple "Cookbook" (CB) scripts demonstrating various use cases of the VBR Calcul
* `CB_011_meltEffects.m` [link](/vbr/examples/CB_011_meltEffects/)
* `CB_012_simplecrust.m` [link](/vbr/examples/CB_012_simplecrust/)
* `CB_013_G_K_inputs.m` [link](/vbr/examples/CB_013_G_K_inputs/)
+* `CB_014_xfit_premelt_extended.m` [link](/vbr/examples/CB_014_xfit_premelt_extended/)
diff --git a/docs/_pages/vbrmethods/anel/xfitpremelt.md b/docs/_pages/vbrmethods/anel/xfitpremelt.md
index 13e6838..4bdfd82 100644
--- a/docs/_pages/vbrmethods/anel/xfitpremelt.md
+++ b/docs/_pages/vbrmethods/anel/xfitpremelt.md
@@ -6,7 +6,9 @@ title: ''
# `xfit_premelt`
-Master curve maxwell scaling using near-solidus parametrization of Yamauchi and Takei (2016), J. Geophys. Res. Solid Earth, [DOI](https://doi.org/10.1002/2016JB013316).
+Master curve maxwell scaling using near-solidus parametrization of Yamauchi and Takei (2016),
+J. Geophys. Res. Solid Earth, [DOI](https://doi.org/10.1002/2016JB013316) with optional extension to include direct melt effects
+of Yamauchi and Takei (2024), J. Geophys. Res. Solid Earth, [DOI](https://doi.org/10.1029/2023JB027738).
## Requires
@@ -21,11 +23,27 @@ VBR.in.SV.sig_MPa % differential stress [MPa]
VBR.in.SV.phi % melt fraction / porosity
VBR.in.SV.rho % density in kg m-3
```
+
+To use the Yamauchi and Takei (2024) scaling that includes direct melt effects, set the
+following flag:
+
+```matlab
+VBR.in.anelastic.xfit_premelt.include_direct_melt_effect = 1;
+```
Additionally, `xfit_premelt` relies on output from the elastic and viscous methods.
-**Required Elastic Methods**: `anharmonic` MUST be in the `VBR.in.elastic.methods_list`. If `anh_poro` is in the methods list then `xfit_premelt` will use the unrelaxed moduli from `anh_poro` (which includes the P,T projection of `anharmonic` plus the poroelastic correction). See the section on [elastic methods](/vbr/vbrmethods/elastic/) for more details.
+**Required Elastic Methods**: `anharmonic` MUST be in the `VBR.in.elastic.methods_list`. Poroelasticity
+is treated differently depending on the value of `include_direct_melt_effect`.
+If `include_direct_melt_effect==0` and `anh_poro` is in the methods list then `xfit_premelt` will use the unrelaxed moduli
+from `anh_poro` (which includes the P,T projection of `anharmonic` plus the poroelastic correction). See the section
+on [elastic methods](/vbr/vbrmethods/elastic/) for more details. If `include_direct_melt_effect==1`, then poroelasticity is incororated within
+J1, following Yamauchi and Takei (2024). The current version of the VBRc uses `include_direct_melt_effect=0` as the default,
+future versions will set this flag to 1 by default.
-**Optional Viscous Methods**: `xfit_premelt` calculates maxwell times using the [viscous xfit_premelt method](/vbr/vbrmethods/visc/xfit_premelt/). If you want to adjust the viscosity calculation used in the maxwell time, you can add `xfit_premelt` to `VBR.in.viscous.methods_list` and adjust the desired parameters. The anelastic calculation will then use the results calculated by the viscous `xfit_premelt` method. This is particularly useful when fitting laboratory measurements of borneol (see [example](/vbr/vbrmethods/anel/xfitpremelt/#example-at-laboratory-conditions) below).
+**Optional Viscous Methods**: `xfit_premelt` calculates maxwell times using the [viscous xfit_premelt method](/vbr/vbrmethods/visc/xfit_premelt/).
+If you want to adjust the viscosity calculation used in the maxwell time, you can add `xfit_premelt` to `VBR.in.viscous.methods_list`
+and adjust the desired parameters. The anelastic calculation will then use the results calculated by the viscous `xfit_premelt` method.
+This is particularly useful when fitting laboratory measurements of borneol (see [example](/vbr/vbrmethods/anel/xfitpremelt/#example-at-laboratory-conditions) below).
## Calling Procedure
@@ -53,6 +71,9 @@ VBR.in.elastic.methods_list={'anharmonic';'anh_poro'};
% set anelastic methods list
VBR.in.anelastic.methods_list={'xfit_premelt'};
+% enable melt effects from Yamauchi and Takei (2024)
+VBR.in.anelastic.xfit_premelt.include_direct_melt_effect = 1;
+
% call VBR_spine
[VBR] = VBR_spine(VBR) ;
```
@@ -95,3 +116,7 @@ The Project script, `Projects/1_LabData/1_Attenuation/FitData_YT16.m` calculates
{:class="img-responsive"}
Data are from figure 10 of Yamauchi and Takei (2016) and are not included in the repository at present.
+
+The Project script, `Projects/1_LabData/1_Attenuation/FitData_YT24.m` reproduces figure 7 from Yamauchi and Takei (2024).
+
+{:class="img-responsive"}
diff --git a/docs/_pages/vbrmethods/anelastic.md b/docs/_pages/vbrmethods/anelastic.md
index bf07170..2d1f1fc 100644
--- a/docs/_pages/vbrmethods/anelastic.md
+++ b/docs/_pages/vbrmethods/anelastic.md
@@ -11,6 +11,9 @@ The anelastic methods displayed by `vbrListMethods` currently include:
* `xfit_mxw` [documentation](/vbr/vbrmethods/anel/xfitmxw/): Master Curve Fit, maxwell scaling
* `xfit_premelt` [documentation](/vbr/vbrmethods/anel/xfitpremelt/): Master Curve Fit, pre-melting maxwell scaling
-A detailed theoretical description of each is provided in the full VBR Manual document (IN PREP) and so here, we only describe the computational aspects useful to the end user.
+A detailed theoretical description of each is provided in the full methods paper ([link](https://doi.org/10.1016/j.pepi.2020.106639)) or
+in the papers cited by the methods and so here, we only describe the computational aspects useful to the end user.
-All of the anelastic methods require results of an elastic calculation, specifically the unrelaxed elastic moduli. If calculated, the anelastic methods will use moduli from `VBR.out.elastic.anh_poro` and default to those from `VBR.out.elastic.anharmonic`.
+All of the anelastic methods require results of an elastic calculation, specifically the unrelaxed elastic moduli.
+If calculated, the anelastic methods may also use moduli from `VBR.out.elastic.anh_poro` and default to those from
+`VBR.out.elastic.anharmonic`.
diff --git a/docs/_pages/vbrmethods/visc/xfit_premelt.md b/docs/_pages/vbrmethods/visc/xfit_premelt.md
index 1afe3f2..7cb8d05 100644
--- a/docs/_pages/vbrmethods/visc/xfit_premelt.md
+++ b/docs/_pages/vbrmethods/visc/xfit_premelt.md
@@ -45,7 +45,11 @@ disp(VBR.in.viscous.xfit_premelt)
### eta_melt_free_method
-The parameter, `VBR.in.viscous.xfit_premelt.eta_melt_free_method`, controls what viscosity method is used to calculate the melt free viscosity. By default, this is set to `xfit_premelt`, in which case it uses the upper mantle fit directly from Yamauchi and Takei (2016). If set to one of the other viscosity mtehods, `HK2003` or `HZK2011`, then the melt free viscosity is calculated using those methods with melt fraction set to 0 and then the near-solidus pre-melting effect is then multiplied on.
+The parameter, `VBR.in.viscous.xfit_premelt.eta_melt_free_method`, controls what viscosity method is used to calculate
+the melt free viscosity. By default, this is set to `xfit_premelt`, in which case it uses the upper mantle fit directly
+from Yamauchi and Takei (2016). If set to one of the other viscosity mtehods, `HK2003` or `HZK2011`, then the melt
+free viscosity is calculated using those methods with melt fraction set to 0 and then the near-solidus pre-melting
+effect is then multiplied on.
## Output
Output is stored in `VBR.out.viscous.xfit_premelt`. Unlike the other viscous methods, `xfit_premelt` only returns a diffusion creep viscosity sub-structure:
diff --git a/docs/assets/images/CBs/CB_004_xfit_premelt.png b/docs/assets/images/CBs/CB_004_xfit_premelt.png
index 5121429..a0d3d4c 100644
Binary files a/docs/assets/images/CBs/CB_004_xfit_premelt.png and b/docs/assets/images/CBs/CB_004_xfit_premelt.png differ
diff --git a/docs/assets/images/CBs/CB_014_xfit_premelt_extended.png b/docs/assets/images/CBs/CB_014_xfit_premelt_extended.png
new file mode 100644
index 0000000..9c162b3
Binary files /dev/null and b/docs/assets/images/CBs/CB_014_xfit_premelt_extended.png differ
diff --git a/docs/assets/images/xfitpremelt_melt_effects.png b/docs/assets/images/xfitpremelt_melt_effects.png
new file mode 100644
index 0000000..059f5c0
Binary files /dev/null and b/docs/assets/images/xfitpremelt_melt_effects.png differ
diff --git a/release_notes.md b/release_notes.md
index 6c7badb..a82e210 100644
--- a/release_notes.md
+++ b/release_notes.md
@@ -1,5 +1,12 @@
# v1.1.3
## New Features
+* updates to xfit_premelt:
+ * add direct melt effects from Yamauchi and Takei, 2024. The `xfit_premelt` method will use the updated parameter values when `VBR.in.anelastic.xfit_premelt.include_direct_melt_effect = 1;` (default is 0, a future VBRc version will change the default to 1).
+ * change default exponential melt factor (the alpha in exp(-alpha*phi) in the viscosity relationship) from 25 to 30.
* add a `VBR_save` function for saving `VBR` structures
-* add framework for handling temporary files in test suite
\ No newline at end of file
+* add framework for handling temporary files in test suite
+* add convenience function, `full_nd`, to create filled arrays
+
+## Bug fix
+* fix for undefined behavior of pre-melt scaling at Tn == 1.0
\ No newline at end of file
diff --git a/vbr/support/vbr_categorical_cmap_array.m b/vbr/support/vbr_categorical_cmap_array.m
new file mode 100644
index 0000000..0b9c193
--- /dev/null
+++ b/vbr/support/vbr_categorical_cmap_array.m
@@ -0,0 +1,39 @@
+function x = vbr_categorical_cmap_array()
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %
+ % x = vbr_categorical_cmap_array()
+ %
+ % return an array of rgb values for a categorical colormap.
+ %
+ % colormap was generated with http://vrl.cs.brown.edu/color , see:
+ % Gramazio et al., 2017, IEEE Transactions on Visualization and Computer
+ % Graphics, "Colorgorical: creating discriminable and preferable color
+ % palettes for information visualization"
+ %
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ x = [get_cat_rgb(88,181,225); ...
+ get_cat_rgb(55,9,35); ...
+ get_cat_rgb(65,210,107); ...
+ get_cat_rgb(177,20,120); ...
+ get_cat_rgb(175,211,90); ...
+ get_cat_rgb(46,14,186); ...
+ get_cat_rgb(120,140,59); ...
+ get_cat_rgb(67,70,171); ...
+ get_cat_rgb(46,236,230); ...
+ get_cat_rgb(0,42,51); ...
+ get_cat_rgb(253,199,204); ...
+ get_cat_rgb(104,60,0); ...
+ get_cat_rgb(192,104,252); ...
+ get_cat_rgb(250,209,57); ...
+ get_cat_rgb(64,96,135); ...
+ get_cat_rgb(35,137,16); ...
+ get_cat_rgb(242,49,252); ...
+ get_cat_rgb(55,78,7); ...
+ get_cat_rgb(161,143,248); ...
+ get_cat_rgb(247,147,2)];
+
+end
+
+function rgb = get_cat_rgb(r, g, b)
+ rgb = [r, g, b] / 255;
+end
\ No newline at end of file
diff --git a/vbr/support/vbr_categorical_color.m b/vbr/support/vbr_categorical_color.m
new file mode 100644
index 0000000..29e09ef
--- /dev/null
+++ b/vbr/support/vbr_categorical_color.m
@@ -0,0 +1,24 @@
+function rgb = vbr_categorical_color(iclr)
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %
+ % rgb = vbr_categorical_color(iclr)
+ %
+ % return a single rgb value from the vbr_categorical_cmap_array.
+ %
+ % Parameters
+ % ----------
+ % iclr
+ % the index to sample from the colormap. Will be wrapped to be within
+ % the bounds of the colormap.
+ %
+ % Output
+ % ------
+ % rgb
+ % 3-element array of floating point rgb values in (0,1) range
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ rgbs = vbr_categorical_cmap_array();
+ ncolors = size(rgbs);
+ ncolors = ncolors(1);
+ iclr = ncolors - mod(iclr, ncolors);
+ rgb = rgbs(iclr, :);
+end
\ No newline at end of file
diff --git a/vbr/testing/test_utilities_full_nd.m b/vbr/testing/test_utilities_full_nd.m
new file mode 100644
index 0000000..0a73dea
--- /dev/null
+++ b/vbr/testing/test_utilities_full_nd.m
@@ -0,0 +1,34 @@
+function TestResult = test_utilities_full_nd()
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% TestResult = test_utilities_full_nd()
+%
+% test full_nd
+%
+% Parameters
+% ----------
+% none
+%
+% Output
+% ------
+% TestResult struct with fields:
+% .passed True if passed, False otherwise.
+% .fail_message Message to display if false
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ x = full_nd(1.0, 2);
+ TestResult.passed = true;
+ TestResult.fail_message = '';
+ if all(size(x)==2) == 0
+ msg = 'full_nd failed full_nd(1.0, 2) check.';
+ TestResult.passed = false;
+ TestResult.fail_message = msg;
+ end
+
+ x = full_nd(1.0, 2, 4, 3);
+ sh = size(x);
+ if sh(1) ~= 2 | sh(2) ~= 4 | sh(3) ~= 3
+ msg = 'full_nd failed full_nd(1.0, 2, 4, 3) check.';
+ TestResult.passed = false;
+ TestResult.fail_message = msg;
+ end
+
+end
\ No newline at end of file
diff --git a/vbr/testing/test_utilities_vbr_categorical_cmap.m b/vbr/testing/test_utilities_vbr_categorical_cmap.m
new file mode 100644
index 0000000..261f964
--- /dev/null
+++ b/vbr/testing/test_utilities_vbr_categorical_cmap.m
@@ -0,0 +1,36 @@
+function TestResult = test_utilities_vbr_categorical_cmap()
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% TestResult = test_utilities_vbr_categorical_cmap()
+%
+% test the categorical colormap
+%
+% Parameters
+% ----------
+% none
+%
+% Output
+% ------
+% TestResult struct with fields:
+% .passed True if passed, False otherwise.
+% .fail_message Message to display if false
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ TestResult.passed = true;
+ TestResult.fail_message = '';
+
+ x = vbr_categorical_cmap_array();
+ if max(x) > 1
+ msg = 'colormap rgb values contain value > 1';
+ TestResult.passed = false;
+ TestResult.fail_message = msg;
+ end
+
+ rgb_1 = vbr_categorical_color(1);
+ rgb_n1 = vbr_categorical_color(numel(x)+1);
+ if all(rgb_1==rgb_n1) == 0
+ msg = 'categorical colormap sample incorrectly wrapped';
+ TestResult.passed = false;
+ TestResult.fail_message = msg;
+ end
+
+end
\ No newline at end of file
diff --git a/vbr/testing/test_vbrcore_007_G_K_TP.m b/vbr/testing/test_vbrcore_007_G_K_TP.m
index 71d73c2..587c2c6 100644
--- a/vbr/testing/test_vbrcore_007_G_K_TP.m
+++ b/vbr/testing/test_vbrcore_007_G_K_TP.m
@@ -70,4 +70,5 @@
VBR.in.elastic.methods_list={'anharmonic';'anh_poro';};
VBR.in.anelastic.methods_list={'eburgers_psp';'andrade_psp';'xfit_mxw'};
VBR.in.viscous.methods_list={'HZK2011'};
+ VBR.in.elastic.quiet = 1;
end
diff --git a/vbr/testing/test_vbrcore_008_premelt_Tn_1.m b/vbr/testing/test_vbrcore_008_premelt_Tn_1.m
new file mode 100644
index 0000000..9c35b0e
--- /dev/null
+++ b/vbr/testing/test_vbrcore_008_premelt_Tn_1.m
@@ -0,0 +1,53 @@
+function TestResult = test_vbrcore_008_premelt_Tn_1()
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% TestResult = test_vbrcore_008_premelt_Tn_1()
+%
+% test for undefined behavior at Tn == 1.0
+%
+% Parameters
+% ----------
+% none
+%
+% Output
+% ------
+% TestResult struct with fields:
+% .passed True if passed, False otherwise.
+% .fail_message Message to display if false
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ TestResult.passed=true;
+ TestResult.fail_message = '';
+
+ %% write method list %%
+ VBR.in.elastic.methods_list={'anharmonic'};
+ VBR.in.anelastic.methods_list={'xfit_premelt'};
+ VBR.in.anelastic.xfit_premelt.include_direct_melt_effect = 0;
+
+ % load anharmonic parameters, adjust Gu_0_ol and derivatives to match YT2016
+ VBR.in.elastic.anharmonic.Gu_0_ol=72.45; %[GPa]
+ VBR.in.elastic.anharmonic.dG_dT = -10.94*1e6; % Pa/C (equivalent ot Pa/K)
+ VBR.in.elastic.anharmonic.dG_dP = 1.987; % GPa / GPa
+
+ %% Define the Thermodynamic State %%
+ VBR.in.SV.T_K=1200:25:1500;
+ VBR.in.SV.T_K=VBR.in.SV.T_K+273;
+ VBR.in.SV.Tsolidus_K = VBR.in.SV.T_K;
+ sz=size(VBR.in.SV.T_K); % temperature [K]
+ VBR.in.SV.P_GPa = full_nd(2.5, sz); % pressure [GPa]
+ VBR.in.SV.dg_um=full_nd(0.01 * 1e6, sz); % grain size [um]
+ VBR.in.SV.rho = full_nd(3300, sz); % density [kg m^-3]
+ VBR.in.SV.sig_MPa = full_nd(.1, sz); % differential stress [MPa]
+ VBR.in.SV.phi = full_nd(0.0, sz);
+ VBR.in.SV.f = 1 ; % 1 Hz
+
+ [VBR] = VBR_spine(VBR);
+
+ if sum(VBR.out.viscous.xfit_premelt.diff.eta == 0.0) > 0
+ msg = 'xfit_premelt.diff.eta contains 0.0 values';
+ TestResult.passed = false;
+ TestResult.fail_message = msg;
+ end
+
+end
+
+
diff --git a/vbr/testing/test_vbrcore_009_premelt_direct.m b/vbr/testing/test_vbrcore_009_premelt_direct.m
new file mode 100644
index 0000000..027e29f
--- /dev/null
+++ b/vbr/testing/test_vbrcore_009_premelt_direct.m
@@ -0,0 +1,55 @@
+function TestResult = test_vbrcore_009_premelt_direct()
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% TestResult = test_vbrcore_009_premelt_direct()
+%
+% test for when include_direct_melt_effect == 1
+%
+% Parameters
+% ----------
+% none
+%
+% Output
+% ------
+% TestResult struct with fields:
+% .passed True if passed, False otherwise.
+% .fail_message Message to display if false
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ TestResult.passed=true;
+ TestResult.fail_message = '';
+
+ %% write method list %%
+ VBR.in.elastic.methods_list={'anharmonic'};
+ VBR.in.anelastic.methods_list={'xfit_premelt'};
+ VBR.in.anelastic.xfit_premelt.include_direct_melt_effect = 1;
+
+ % load anharmonic parameters, adjust Gu_0_ol and derivatives to match YT2016
+ VBR.in.elastic.anharmonic.Gu_0_ol=72.45; %[GPa]
+ VBR.in.elastic.anharmonic.dG_dT = -10.94*1e6; % Pa/C (equivalent ot Pa/K)
+ VBR.in.elastic.anharmonic.dG_dP = 1.987; % GPa / GPa
+
+ %% Define the Thermodynamic State %%
+ VBR.in.SV.phi = logspace(-8,log10(0.05), 10);
+
+ sz=size(VBR.in.SV.phi); % temperature [K]
+ VBR.in.SV.T_K = full_nd(1350+273, sz);
+ VBR.in.SV.Tsolidus_K = full_nd(1300+273., sz);
+ VBR.in.SV.P_GPa = full_nd(2.5, sz); % pressure [GPa]
+ VBR.in.SV.dg_um=full_nd(0.01 * 1e6, sz); % grain size [um]
+ VBR.in.SV.rho = full_nd(3300, sz); % density [kg m^-3]
+ VBR.in.SV.sig_MPa = full_nd(.1, sz); % differential stress [MPa]
+
+ VBR.in.SV.f = 10.; % 1 Hz
+
+ [VBR] = VBR_spine(VBR);
+
+ Qinv = VBR.out.anelastic.xfit_premelt.Qinv;
+ phi = VBR.in.SV.phi;
+ dQdphi = (Qinv(2:end) - Qinv(1:end-1))./ (phi(2:end) - phi(1:end-1));
+ if sum(dQdphi<0) > 0
+ TestResult.passed=false;
+ TestResult.fail_message = 'Increasing melt fraction should increase Qinv.';
+ end
+end
+
+
diff --git a/vbr/vbrCore/functions/Q_xfit_premelt.m b/vbr/vbrCore/functions/Q_xfit_premelt.m
index 62f637b..c549e2e 100644
--- a/vbr/vbrCore/functions/Q_xfit_premelt.m
+++ b/vbr/vbrCore/functions/Q_xfit_premelt.m
@@ -16,28 +16,50 @@
disp('To use Q_xfit_premelt, you must provide VBR.in.SV.Tsolidus_K')
end
+ params=VBR.in.anelastic.xfit_premelt;
+
if has_solidus
% state variables
- if isfield(VBR.in.elastic,'anh_poro')
- Gu_in = VBR.out.elastic.anh_poro.Gu;
- elseif isfield(VBR.in.elastic,'anharmonic')
- Gu_in = VBR.out.elastic.anharmonic.Gu;
+ if params.include_direct_melt_effect == 1
+ % in YT2024, all poroealstic effects are applied internally to J1, so
+ % always use anharmonic as unrelaxed
+ Gu_in = VBR.out.elastic.anharmonic.Gu;
+ else
+ if isfield(VBR.in.elastic,'anh_poro')
+ Gu_in = VBR.out.elastic.anh_poro.Gu;
+ elseif isfield(VBR.in.elastic,'anharmonic')
+ Gu_in = VBR.out.elastic.anharmonic.Gu;
+ end
end
+
Ju_in = 1./Gu_in ;
rho = VBR.in.SV.rho ;
phi = VBR.in.SV.phi ;
Tn=VBR.in.SV.T_K./VBR.in.SV.Tsolidus_K ; % solidus-normalized temperature
- params=VBR.in.anelastic.xfit_premelt;
% maxwell time
[tau_m,VBR]=MaxwellTimes(VBR,Gu_in);
% calculate the Tn-dependent coefficients, A_p and sig_p
[A_p,sig_p]=calcApSigp(Tn,phi,params);
+ if params.include_direct_melt_effect == 1
+ Beta_B = params.Beta_B * phi;
+ else
+ Beta_B = 0.0;
+ end
+
+ % poroelastic J1 effect if applicable
+ if params.include_direct_melt_effect == 1
+ % poroelastic effect added to J1, Delta_poro
+ poro_elastic_factor = params.poro_Lambda * phi;
+ else
+ % no poroelastic effects outside of incoming unrelaxed modulus
+ poro_elastic_factor = 0.0;
+ end
% set other constants
alpha_B=params.alpha_B;
- A_B=params.A_B;
+ A_B_plus_Beta_B= params.A_B + Beta_B;
tau_pp=params.tau_pp;
% set up frequency dependence
@@ -48,20 +70,21 @@
n_freq = numel(VBR.in.SV.f);
sz = size(Gu_in);
J1 = proc_add_freq_indeces(zeros(sz),n_freq);
- J2 = J1; V = J1;
+ J2 = J1; V = J1; f_norm_glob=J1;
n_SVs=numel(Tn); % total elements in state variables
% loop over frequencies, calculate J1,J2
for i = 1:n_freq
+
sv_i0=(i-1)*n_SVs + 1; % starting linear index of this freq
sv_i1=sv_i0+n_SVs-1; % ending linear index of this freq
-
+ f_norm_glob(sv_i0:sv_i1)=tau_m*VBR.in.SV.f(i);
p_p=period_vec(i)./(2*pi*tau_m);
% tau_eta^S= tau_s / (2 pi tau_m);, tau_s = seismic wave period, tau_m = ss maxwell time
- ABppa=A_B*(p_p.^alpha_B);
+ ABppa=A_B_plus_Beta_B .* (p_p.^alpha_B);
lntaupp=log(tau_pp./p_p);
- J1(sv_i0:sv_i1)=Ju_in .* (1 + ABppa/alpha_B+ ...
+ J1(sv_i0:sv_i1)=Ju_in .* (1 + poro_elastic_factor + ABppa/alpha_B+ ...
pifac*A_p.*sig_p.*(1-erf(lntaupp./(sqrt(2).*sig_p))));
J2(sv_i0:sv_i1)=Ju_in*pi/2.*(ABppa+A_p.*(exp(-(lntaupp.^2)./(2*sig_p.^2)))) + ...
Ju_in.*p_p;
@@ -80,6 +103,7 @@
VBRout.Qinv = J2./J1.*(J2_J1_frac.^-1);
VBRout.Q = 1./VBRout.Qinv;
VBRout.M=1./sqrt(J1.^2+J2.^2);
+ VBRout.f_norm = f_norm_glob;
% calculate mean velocity along frequency dimension
VBRout.Vave = Q_aveVoverf(VBRout.V,VBR.in.SV.f);
@@ -87,6 +111,7 @@
VBRout.units = Q_method_units();
VBRout.units.M1 = 'Pa';
VBRout.units.M2 = 'Pa';
+ VBRout.units.f_norm = '';
% store the output structure
VBR.out.anelastic.xfit_premelt=VBRout;
@@ -101,11 +126,15 @@
Ap_Tn_pts=params.Ap_Tn_pts;
sig_p_Tn_pts=params.sig_p_Tn_pts;
- Beta=params.Beta; %
+ if params.include_direct_melt_effect == 1
+ Beta_p = params.Beta; % this depends on melt fraction in YT2024
+ else
+ Beta_p = 0; % no direct melt dependence
+ end
A_p=zeros(size(Tn));
- A_p(Tn>=Ap_Tn_pts(3))=params.Ap_fac_3+Beta*phi(Tn>=Ap_Tn_pts(3));
- A_p(Tn= Ap_Tn_pts(3))=params.Ap_fac_3+Beta_p*phi(Tn >= Ap_Tn_pts(3));
+ A_p(Tn < Ap_Tn_pts(3))=params.Ap_fac_3;
+ A_p(Tn < Ap_Tn_pts(2))=params.Ap_fac_1 + params.Ap_fac_2*(Tn(Tn < Ap_Tn_pts(2))-Ap_Tn_pts(1));
A_p(Tn < Ap_Tn_pts(1))=params.Ap_fac_1;
sig_p=zeros(size(Tn));
diff --git a/vbr/vbrCore/functions/el_ModUnrlx_dTdP_f.m b/vbr/vbrCore/functions/el_ModUnrlx_dTdP_f.m
index fdb7a2d..5d8f338 100644
--- a/vbr/vbrCore/functions/el_ModUnrlx_dTdP_f.m
+++ b/vbr/vbrCore/functions/el_ModUnrlx_dTdP_f.m
@@ -37,9 +37,19 @@
Gu_TP = VBR.in.elastic.Gu_TP; % Pa
% calculate bulk modulus
-disp(['Unrelaxed shear modulus was provided without an unrelaxed bulk modulus.', ...
- ' Calculating bulk modulus assuming nu=', num2str(nu), ...
- '. Set VBR.in.elastic.Ku_TP to specify a value.']);
+ if isfield(VBR.in.elastic,'quiet')
+ quiet_mode = VBR.in.elastic.quiet;
+ else
+ quiet_mode = 0;
+ end
+
+ if quiet_mode == 0
+ msg = ['Unrelaxed shear modulus was provided without an unrelaxed ', ...
+ 'bulk modulus. Calculating bulk modulus assuming nu=', num2str(nu), ...
+ '. Set VBR.in.elastic.Ku_TP to specify a value or set ', ...
+ 'VBR.in.elastic.quiet = 1 to silence this warning.'];
+ disp(msg);
+ end
Ku_TP = calc_Ku(Gu_TP,nu);
else
dT = (VBR.in.SV.T_K-T_K_ref);
diff --git a/vbr/vbrCore/functions/full_nd.m b/vbr/vbrCore/functions/full_nd.m
new file mode 100644
index 0000000..4abb0bb
--- /dev/null
+++ b/vbr/vbrCore/functions/full_nd.m
@@ -0,0 +1,23 @@
+function X = full_nd(fill_val, varargin)
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %
+ % X = full_nd(fill_val, N)
+ %
+ % returns an N-D array filled with a constant. After fill_val, all arguments
+ % are forwarded to ones().
+ %
+ % Parameters:
+ % ----------
+ % fill_val
+ % the number to fill the array (or matrix) with
+ %
+ % remaining arguments are forwarded to ones().
+ %
+ %
+ % Output:
+ % ------
+ % matrix
+ %
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ X = ones(varargin{:}) * fill_val;
+end
\ No newline at end of file
diff --git a/vbr/vbrCore/functions/visc_calc_xfit_premelt.m b/vbr/vbrCore/functions/visc_calc_xfit_premelt.m
index 35cd6cc..1d8ddcc 100644
--- a/vbr/vbrCore/functions/visc_calc_xfit_premelt.m
+++ b/vbr/vbrCore/functions/visc_calc_xfit_premelt.m
@@ -130,7 +130,7 @@
T_eta=params.T_eta;
gamma=params.gamma;
lambda=params.alpha; % rename to be consistent with YT2016 nomenclature
-
+ B = params.B;
A_n=zeros(size(Tn));
A_n(Tn= T_eta) & (Tn < 1);
A_n(msk)=exp(-(Tn(msk)-T_eta)./(Tn(msk)-Tn(msk)*T_eta)*log(gamma));
- msk=(Tn > 1);
- A_n(msk)=exp(-lambda*phi(msk))/gamma;
+ msk=(Tn >= 1);
+ A_n(msk)=exp(-lambda*phi(msk))/gamma/B;
end
diff --git a/vbr/vbrCore/params/Params_Anelastic.m b/vbr/vbrCore/params/Params_Anelastic.m
index 62d7e78..dd41324 100644
--- a/vbr/vbrCore/params/Params_Anelastic.m
+++ b/vbr/vbrCore/params/Params_Anelastic.m
@@ -102,7 +102,8 @@
if strcmp(method,'xfit_premelt')
% xfit_premelt parameters
- params.citations={'Yamauchi and Takei, 2016, J. Geophys. Res. Solid Earth, https://doi.org/10.1002/2016JB013316'};
+ params.citations={'Yamauchi and Takei, 2016, J. Geophys. Res. Solid Earth, https://doi.org/10.1002/2016JB013316';
+ 'Yamauchi and Takei, 2024, J. Geophys. Res. Solid Earth, https://doi.org/10.1029/2023JB027738'};
params.func_name='Q_xfit_premelt'; % the name of the matlab function
% high temp background spectrum
@@ -110,8 +111,7 @@
params.A_B=0.664; % high temp background dissipation strength
% pre-melting dissipation peak settings:
- params.tau_pp=6*1e-5; % peak center
- params.Beta=0; %
+ params.tau_pp=6*1e-5; % peak center, table 4 of YT16, paragraph before eq 10
params.Ap_fac_1=0.01;
params.Ap_fac_2=0.4;
params.Ap_fac_3=0.03;
@@ -120,6 +120,20 @@
params.sig_p_fac_3=7;
params.Ap_Tn_pts=[0.91,0.96,1]; % Tn cuttoff points
params.sig_p_Tn_pts=[0.92,1]; % Tn cuttoff points
+
+ % note: sig_p_fac2 and Ap_fac_2 above are combinations of the other constants
+ % and appear in the original papers as:
+ % sig_p_fac2 = (sig_p_fac_3-sig_p_fac_1) / (sig_p_Tn_pts(2) - sig_p_Tn_pts(1))
+ % Ap_fac_2 = (Ap_fac_3-Ap_fac_1) / (Ap_Tn_pts(3) - Ap_Tn_pts(1))
+
+ % melt effects. The following beta values are set to 0.0 within Q_xfit_premelt
+ % if include_direct_melt_effect = 0, corresponding to YT2016. If set to 1,
+ % the scaling will follow YT2024. Additionally, include_direct_melt_effect=1
+ % will trigger different poro-elastic behavior.
+ params.include_direct_melt_effect = 0; % set to 1 to include YT2024 melt effect
+ params.Beta=1.38; % this is determined in YT2024, named Beta_P in YT2024 eq 5
+ params.Beta_B=6.94; % YT2024 only
+ params.poro_Lambda = 4.0; % Table 6 YT2024,
params.description='pre-melting scaling';
end
diff --git a/vbr/vbrCore/params/Params_Viscous.m b/vbr/vbrCore/params/Params_Viscous.m
index fd5dc27..d33c266 100644
--- a/vbr/vbrCore/params/Params_Viscous.m
+++ b/vbr/vbrCore/params/Params_Viscous.m
@@ -47,9 +47,10 @@
params.citations={'Yamauchi and Takei, 2016, J. Geophys. Res. Solid Earth, https://doi.org/10.1002/2016JB013316'};
params.description='Steady state flow law for pre-melting viscosity drop, flow law parameters from fitting upper mantle.';
% near-solidus and melt effects
- params.alpha=25; % taken from diff. creep value of HZK2011. YT2016 call this lambda.
+ params.alpha=30; % taken from diff. creep value of HZK2011. YT2016 call this lambda.
params.T_eta=0.94;
params.gamma=5;
+ params.B = 1.0;
% method to use for melt-free diff. creep viscosity
params.eta_melt_free_method='xfit_premelt';