MATLAB scripts for "On the influence of mineral dust on glacial albedo at Nevado Huascarán"

On the influence of mineral dust on glacial albedo at Nevado Huascarán Austin M. Weber | Byrd Polar and Climate Research Center, Columbus, Ohio, USA Note on dependencies: This script requires functions from Toolboxes and File Exchange submissions that are not a part of the base MATLAB language. You will need: • Statistics and Machine Learning Toolbox • Curve Fitting Toolbox • Climate Data Toolbox (Greene et al., 2019) • Categorical Colors (Weber, 2025) Define the path to the data datafile = 'dust-data-aggregated-scnicar.xlsx'; Load the data colB_5y = readDustTable(datafile,'ColB_5y'); % [ug dust / g snow] summitB_5y = readDustTable(datafile,'SummitB_5y'); % [ug dust / g snow] colB_SNICAR = readDustTable(datafile,'ColB_SNICAR'); % [Broadband albedo] summitB_SNICAR = readDustTable(datafile,'SummitB_SNICAR'); % [Broadband albedo] Compare total dust records age5y = colB_5y.StartYear; age5y_labels = {'1960-4','1965-9','1970-4','1975-9','1980-4','1985-9',... '1990-4','1995-9','2000-4','2005-9','2010-4','2015-9'}; colB_5y_totalDust = sum(colB_5y{:,3:end},2); summitB_5y_totalDust = sum(summitB_5y{:,3:end},2); lw=3; tiledlayout(2,2,'TileSpacing','compact') nexttile([1 2]) s1=stairs(age5y,[colB_5y_totalDust summitB_5y_totalDust],... 'LineWidth',lw,... 'Color','k'); s1(2).LineStyle=':'; ylim([0 5]) ylabel('Total Dust (µg/g)') xticks(1960:5:2015) xlim([1959 2016]) xticklabels(age5y_labels) lg=legend('Col','Summit','Location','northwest'); lg.IconColumnWidth=15; 1 lg.Box='off'; updateFont updateGrid nexttile() stairs(age5y,colB_5y{:,3:end},... 'LineWidth',lw) xticks(1960:5:2015) xlim([1959 2016]) xticklabels(age5y_labels) ylabel('Col (µg/g)') lg2=legend('0.63-1.00 µm',... '1.00-2.52 µm',... '2.52-5.04 µm',... '5.04-10.07 µm',... '10.07-16.00 µm',... 'Location','northwest'); lg2.IconColumnWidth=12; lg2.Box='off'; updateFont updateGrid updateColorOrder nexttile stairs(age5y,summitB_5y{:,3:end},':',... 'LineWidth',lw) xticks(1960:5:2015) xlim([1959 2016]) xticklabels(age5y_labels) ylabel('Summit (µg/g)') updateFont updateGrid updateColorOrder 2 Compare former (1960-1989) and latter (1990-2019) periods col = repmat({'Col'},[12 1]); summit = repmat({'Summit'},[12 1]); former = repmat(1960,[6 1]); latter = repmat(1990,[6 1]); formerLatterTable = array2table([col; summit],'VariableNames',{'Site'}); formerLatterTable.TotalDust = [colB_5y_totalDust; summitB_5y_totalDust]; formerLatterTable.Period = [former; latter; former; latter]; %%% %%% Visualize %%% figure groupedSpacedBoxchart(formerLatterTable,'Site','TotalDust','Period'); colororder(cat_colors('colorblind')) yticks(0.5:0.5:5),yticklabels({'0.5','1.0','1.5','2.0','2.5','3.0','3.5','4.0','4.5' ,'5.0'}) ylabel('Concentration (µg/g)') lg=legend('1960-1989','1990-2019'); lg.Box='off'; lg.IconColumnWidth=14; updateFont updateGrid box on axis square 3 %%% %%% Evaluate statistical difference scores %%% % Test whether the means are statistically different h = ttest2(colB_5y_totalDust(1:6),colB_5y_totalDust(7:12)); if h == 1 fprintf('The null hypothesis is rejected. The COL sample means are statistically different.') else fprintf('The null hypothesis cannot be rejected. The COL sample means are not statistically different.') end The null hypothesis is rejected. The COL sample means are statistically different. % Test whether the medians are statistically different p = ranksum(colB_5y_totalDust(1:6),colB_5y_totalDust(7:12)); if p < 0.05 fprintf('The COL sample medians are statistically different.') else fprintf('The COL sample medians are not statistically different.') end The COL sample medians are statistically different. % Test whether the means are statistically different h = ttest2(summitB_5y_totalDust(1:6),summitB_5y_totalDust(7:12)); 4 if h == 1 fprintf('The null hypothesis is rejected. The SUMMIT sample means are statistically different.') else fprintf('The null hypothesis cannot be rejected. The SUMMIT sample means are not statistically different.') end The null hypothesis is rejected. The SUMMIT sample means are statistically different. % Test whether the medians are statistically different p = ranksum(summitB_5y_totalDust(1:6),summitB_5y_totalDust(7:12)); if p < 0.05 fprintf('The SUMMIT sample medians are statistically different.') else fprintf('The SUMMIT sample medians are not statistically different.') end The SUMMIT sample medians are statistically different. SNICAR-ADv3 Modeled Albedos lw=2.5; xtickvals = 1960:5:2015; years = colB_SNICAR.Sample; close figure % Synthetic data for organizing the legend--these datapoints will not be shown % on the chart synthX = [1 2 3]; synthY = [-1 -0.5 0]; plot(synthX,synthY,'k-','LineWidth',lw); hold on plot(synthX,synthY,'k:','LineWidth',lw); plot(synthX,synthY,'k','Marker','o','MarkerSize',8,'MarkerFaceColor','w','LineStyle' ,'none'); plot(synthX,synthY,'k','Marker','diamond','MarkerSize',8,'MarkerFaceColor','w','Line Style','none'); hold on % Control ylim([0.714 0.734]) yline(0.732956,'-r','linewidth',lw) % Col data col = plot(xtickvals,colB_SNICAR{:,2:end},... 'Color','k',... 'LineWidth',lw); col(1).Marker='o'; col(1).MarkerSize=8; col(1).MarkerFaceColor='w'; col(2).Marker='diamond'; col(2).MarkerSize=8; col(2).MarkerFaceColor='w'; xlim([1958 2017]) xticks(1960:5:2015) xticklabels(years) 5 ylabel('Broadband Albedo') hold on % Summit data summit = plot(xtickvals,summitB_SNICAR{:,2:end},':',... 'Color','k',... 'LineWidth',lw); summit(1).Marker='o'; summit(1).MarkerSize=8; summit(1).MarkerFaceColor='w'; summit(2).Marker='diamond'; summit(2).MarkerSize=8; summit(2).MarkerFaceColor='w'; %ylim([0.8259 0.8306]); yticks(0.714 : 0.002 : 0.734) yticklabels({'0.714','0.716','0.718','0.720','0.722','0.724',... '0.726','0.728','0.730','0.732','0.734'}) lg = legend('Col','Summit','Dust','Dust + BC','Control',... 'Location','northeast'); lg.IconColumnWidth=14; lg.Box='off'; updateFont updateGrid Compare former (1960-1989) and latter (1990-2019) periods col = repmat({'Col'},[12 1]); summit = repmat({'Summit'},[12 1]); former = repmat(1960,[6 1]); latter = repmat(1990,[6 1]); formerLatterAlbedoTable = array2table([col; summit],'VariableNames',{'Site'}); 6 formerLatterAlbedoTable.TotalDust = [colB_SNICAR.Albedo; summitB_SNICAR.Albedo]; formerLatterAlbedoTable.Period = [former; latter; former; latter]; %%% %%% Visualize %%% figure groupedSpacedBoxchart(formerLatterAlbedoTable,'Site','TotalDust','Period'); colororder(cat_colors('colorblind')) ylim([0.714 0.728]) yticks(0.714:0.002:0.728), yticklabels({'0.714','0.716','0.718',... '0.720','0.722','0.724','0.726','0.728'}) ylabel('Broadband Albedo') lg=legend('1960-1989','1990-2019',... 'Location','northwest'); lg.Box='off'; lg.IconColumnWidth=14; updateFont updateGrid box on axis square %%% %%% Evaluate statistical difference scores %%% % Test whether the means are statistically different 7 h = ttest2(colB_SNICAR.Albedo(1:6),colB_SNICAR.Albedo(7:12)); if h == 1 fprintf('The null hypothesis is rejected. The COL sample means are statistically different.') else fprintf('The null hypothesis cannot be rejected. The COL sample means are not statistically different.') end The null hypothesis is rejected. The COL sample means are statistically different. % Test whether the medians are statistically different p = ranksum(colB_SNICAR.Albedo(1:6),colB_SNICAR.Albedo(7:12)); if p < 0.05 fprintf('The COL sample medians are statistically different.') else fprintf('The COL sample medians are not statistically different.') end The COL sample medians are statistically different. % Test whether the means are statistically different h = ttest2(summitB_SNICAR.Albedo(1:6),summitB_SNICAR.Albedo(7:12)); if h == 1 fprintf('The null hypothesis is rejected. The SUMMIT sample means are statistically different.') else fprintf('The null hypothesis cannot be rejected. The SUMMIT sample means are not statistically different.') end The null hypothesis is rejected. The SUMMIT sample means are statistically different. % Test whether the medians are statistically different p = ranksum(summitB_SNICAR.Albedo(1:6),summitB_SNICAR.Albedo(7:12)); if p < 0.05 fprintf('The SUMMIT sample medians are statistically different.') else fprintf('The SUMMIT sample medians are not statistically different.') end The SUMMIT sample medians are statistically different. Logarithmic fitting %%% %%% COL %%% % Get data since 1960 yrs = (1960:5:2015)'; col_albedos = colB_SNICAR.Albedo; % Fit model 8 [f_col,s_col] = fit(yrs,col_albedos,"log10") f_col = Linear model Log10: f_col(x) = a*log10(x) + b Coefficients (with 95% confidence bounds): a = -0.4948 (-0.7263, -0.2633) b = 2.352 (1.588, 3.115) s_col = struct with fields: sse: 0.00 rsquare: 0.69 dfe: 10.00 adjrsquare: 0.66 rmse: 0.00 col_logMdl = @(x) -0.4948*log10(x)+2.352; % Generate the logarithmic curve and confidence intervals x_fit = 1955:1:2105; [col_ci, col_fit] = predint(f_col, x_fit, 0.95); %%% %%% SUMMIT %%% summit_albedos = summitB_SNICAR.Albedo; [f_summit,s_summit] = fit(yrs,summit_albedos,"log10") f_summit = Linear model Log10: f_summit(x) = a*log10(x) + b Coefficients (with 95% confidence bounds): a = -0.4491 (-0.7111, -0.1871) b = 2.204 (1.34, 3.068) s_summit = struct with fields: sse: 0.00 rsquare: 0.59 dfe: 10.00 adjrsquare: 0.55 rmse: 0.00 summit_logMdl = @(x) -0.4491*log10(x)+2.204; [summit_ci, summit_fit] = predint(f_summit, x_fit, 0.95); %%% %%% VISUALIZE %%% scalingFactor = 1.2; close figure('units','inches','position',[1 1 7/scalingFactor 4/scalingFactor]) % Col boundedline(x_fit,col_fit,[col_fit-col_ci(:,1), col_ci(:,2)-col_fit],'- k','alpha','linewidth',3) ylabel('Broadband Albedo') hold on % Summit boundedline(x_fit,summit_fit,[summit_fit-summit_ci(:,1), summit_ci(:,2)- summit_fit],':k','alpha','linewidth',3) 9