diff --git a/Projects/vbr_core_examples/CB_002_2D_HalfSpaceCooling.m b/Projects/vbr_core_examples/CB_002_2D_HalfSpaceCooling.m
index 479e452..51fb6f5 100644
--- a/Projects/vbr_core_examples/CB_002_2D_HalfSpaceCooling.m
+++ b/Projects/vbr_core_examples/CB_002_2D_HalfSpaceCooling.m
@@ -85,7 +85,7 @@
      title(['V_s [km/s] andrade_psp at ',num2str(VBR.in.SV.f(i_f)),' Hz'])
      colorbar()
   end
-
+  saveas(gcf,'./figures/CB_002_2D_HalfSpaceCooling.png')
   % contour percent difference in shear wave velo between two anelastic methods
   % at different frequencies
   dV=abs(VBR.out.anelastic.andrade_psp.V-VBR.out.anelastic.xfit_mxw.V);
diff --git a/Projects/vbr_core_examples/CB_004_xfit_premelt.m b/Projects/vbr_core_examples/CB_004_xfit_premelt.m
index 09eac66..c4d1951 100644
--- a/Projects/vbr_core_examples/CB_004_xfit_premelt.m
+++ b/Projects/vbr_core_examples/CB_004_xfit_premelt.m
@@ -60,3 +60,4 @@
   subplot(1,3,3)
   semilogx(1./VBR.in.SV.f,1e-3*squeeze(VBR.out.anelastic.xfit_premelt.V(1,:,:)));
   ylabel('V_s [km/s]'); xlabel('period [s]')
+  saveas(gcf,'./figures/CB_004_xfit_premelt.png')
diff --git a/Projects/vbr_core_examples/CB_006_viscosity.m b/Projects/vbr_core_examples/CB_006_viscosity.m
index 9a00f1f..339d701 100644
--- a/Projects/vbr_core_examples/CB_006_viscosity.m
+++ b/Projects/vbr_core_examples/CB_006_viscosity.m
@@ -68,3 +68,4 @@
 
   subplot(2,2,4)
   box on; xlabel('log10 \sigma [MPa]'); ylabel('log10 effective viscosity [Pa s]')
+saveas(gcf,'./figures/CB_006_viscosity.png')
diff --git a/Projects/vbr_core_examples/CB_007_SLB2005.m b/Projects/vbr_core_examples/CB_007_SLB2005.m
index 0dd51e3..c0b6d37 100644
--- a/Projects/vbr_core_examples/CB_007_SLB2005.m
+++ b/Projects/vbr_core_examples/CB_007_SLB2005.m
@@ -43,3 +43,4 @@
   plot(VBR.out.elastic.SLB2005.Vs,z/1e3)
   set(gca,'Ydir','reverse')
   xlabel('Vs [km/s]'); ylabel('z [km]')
+  saveas(gcf,'./figures/CB_007_SLB2005.png')
diff --git a/Projects/vbr_core_examples/CB_008_anharmonic_Guo.m b/Projects/vbr_core_examples/CB_008_anharmonic_Guo.m
index d4487cb..833185c 100644
--- a/Projects/vbr_core_examples/CB_008_anharmonic_Guo.m
+++ b/Projects/vbr_core_examples/CB_008_anharmonic_Guo.m
@@ -69,3 +69,5 @@
     ylim([0,150])
     set(gca,'ydir','reverse')
   end
+
+  saveas(gcf,'./figures/CB_008_anharmonic_Gu0.png')
diff --git a/Projects/vbr_core_examples/CB_009_anhporo.m b/Projects/vbr_core_examples/CB_009_anhporo.m
index dc8ef26..513460d 100644
--- a/Projects/vbr_core_examples/CB_009_anhporo.m
+++ b/Projects/vbr_core_examples/CB_009_anhporo.m
@@ -34,3 +34,5 @@
   semilogx(VBR.in.SV.phi,VBR.out.elastic.anh_poro.Vsu/1e3,'k','linewidth',1.5)
   xlabel('\phi'); ylabel('Vsu(P,T,\phi) [km/s]')
   set(gca,'linewidth',1.5)
+
+  saveas(gcf,'./figures/CB_009_anhporo.png')
diff --git a/Projects/vbr_core_examples/CB_010_depthProfiles.m b/Projects/vbr_core_examples/CB_010_depthProfiles.m
index 9d30aca..8e10968 100644
--- a/Projects/vbr_core_examples/CB_010_depthProfiles.m
+++ b/Projects/vbr_core_examples/CB_010_depthProfiles.m
@@ -50,7 +50,7 @@
 
   %% Build figures %%
 
-    figure()
+    figure('PaperPosition',[0,0,8,4],'PaperPositionMode','manual')
     ax1=subplot(1,5,1);
     plot(HF.T_C,HF.z_km)
     hold on
@@ -93,6 +93,7 @@
       xlim([0,8])
     end
 
+    saveas(gcf,'./figures/CB_00210_depthProfiles.png')
 end
 
 function HF = HalfspaceModel(age_Myrs)
diff --git a/Projects/vbr_core_examples/CB_013_G_K_inputs.m b/Projects/vbr_core_examples/CB_013_G_K_inputs.m
new file mode 100644
index 0000000..9bb81fd
--- /dev/null
+++ b/Projects/vbr_core_examples/CB_013_G_K_inputs.m
@@ -0,0 +1,35 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% CB_013_G_K_inputs.m
+%
+%  Specify unrelaxed shear and bulk moduli.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% put the VBRc in the path %
+clear
+path_to_top_level_vbr='../../';
+addpath(path_to_top_level_vbr)
+vbr_init
+
+% specify state variables as usual
+VBR = struct();
+VBR.in.SV.T_K = linspace(800, 1000, 4);
+sz_T = size(VBR.in.SV.T_K);
+VBR.in.SV.P_GPa = linspace(2, 3, 4);
+VBR.in.SV.rho = 3300 * ones(sz_T);
+VBR.in.SV.phi = 0.01 * ones(sz_T);
+VBR.in.SV.sig_MPa = 1 * ones(sz_T);
+VBR.in.SV.dg_um = 1e4 * ones(sz_T);
+VBR.in.SV.f = [0.01, 0.1];
+
+% specify methods as usual
+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'};
+
+% also specify the unrelaxed moduli at elevated temperature, pressure.
+% you could, instead, load these from a file!
+VBR.in.elastic.Gu_TP = linspace(50, 60, 4) * 1e9; % shear modulus
+VBR.in.elastic.Ku_TP = VBR.in.elastic.Gu_TP * 1.5; % bulk modulus
+
+% call the VBRc
+VBR = VBR_spine(VBR);
diff --git a/docs/_pages/examples/CB_002_2D_HalfSpaceCooling.md b/docs/_pages/examples/CB_002_2D_HalfSpaceCooling.md
index 8cfd4ce..8e90741 100644
--- a/docs/_pages/examples/CB_002_2D_HalfSpaceCooling.md
+++ b/docs/_pages/examples/CB_002_2D_HalfSpaceCooling.md
@@ -4,6 +4,9 @@ title: ""
 ---
 
 # CB_002_2D_HalfSpaceCooling.m
+## output figures
+
+!['CB_002_2D_HalfSpaceCooling'](/vbr/assets/images/CBs/CB_002_2D_HalfSpaceCooling.png){:class="img-responsive"}
 ## contents
 ```matlab
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -93,7 +96,7 @@ title: ""
      title(['V_s [km/s] andrade_psp at ',num2str(VBR.in.SV.f(i_f)),' Hz'])
      colorbar()
   end
-
+  saveas(gcf,'./figures/CB_002_2D_HalfSpaceCooling.png')
   % contour percent difference in shear wave velo between two anelastic methods
   % at different frequencies
   dV=abs(VBR.out.anelastic.andrade_psp.V-VBR.out.anelastic.xfit_mxw.V);
diff --git a/docs/_pages/examples/CB_004_xfit_premelt.md b/docs/_pages/examples/CB_004_xfit_premelt.md
index ef38e96..d3826b8 100644
--- a/docs/_pages/examples/CB_004_xfit_premelt.md
+++ b/docs/_pages/examples/CB_004_xfit_premelt.md
@@ -4,6 +4,9 @@ title: ""
 ---
 
 # CB_004_xfit_premelt.m
+## output figures
+
+!['CB_004_xfit_premelt'](/vbr/assets/images/CBs/CB_004_xfit_premelt.png){:class="img-responsive"}
 ## contents
 ```matlab
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -68,4 +71,5 @@ title: ""
   subplot(1,3,3)
   semilogx(1./VBR.in.SV.f,1e-3*squeeze(VBR.out.anelastic.xfit_premelt.V(1,:,:)));
   ylabel('V_s [km/s]'); xlabel('period [s]')
+  saveas(gcf,'./figures/CB_004_xfit_premelt.png')
 ```
diff --git a/docs/_pages/examples/CB_006_viscosity.md b/docs/_pages/examples/CB_006_viscosity.md
index 7aa87e6..18a8208 100644
--- a/docs/_pages/examples/CB_006_viscosity.md
+++ b/docs/_pages/examples/CB_006_viscosity.md
@@ -4,6 +4,9 @@ title: ""
 ---
 
 # CB_006_viscosity.m
+## output figures
+
+!['CB_006_viscosity'](/vbr/assets/images/CBs/CB_006_viscosity.png){:class="img-responsive"}
 ## contents
 ```matlab
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -76,4 +79,5 @@ title: ""
 
   subplot(2,2,4)
   box on; xlabel('log10 \sigma [MPa]'); ylabel('log10 effective viscosity [Pa s]')
+saveas(gcf,'./figures/CB_006_viscosity.png')
 ```
diff --git a/docs/_pages/examples/CB_007_SLB2005.md b/docs/_pages/examples/CB_007_SLB2005.md
index ba0fcc2..4e41f45 100644
--- a/docs/_pages/examples/CB_007_SLB2005.md
+++ b/docs/_pages/examples/CB_007_SLB2005.md
@@ -4,6 +4,9 @@ title: ""
 ---
 
 # CB_007_SLB2005.m
+## output figures
+
+!['CB_007_SLB2005'](/vbr/assets/images/CBs/CB_007_SLB2005.png){:class="img-responsive"}
 ## contents
 ```matlab
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -51,4 +54,5 @@ title: ""
   plot(VBR.out.elastic.SLB2005.Vs,z/1e3)
   set(gca,'Ydir','reverse')
   xlabel('Vs [km/s]'); ylabel('z [km]')
+  saveas(gcf,'./figures/CB_007_SLB2005.png')
 ```
diff --git a/docs/_pages/examples/CB_008_anharmonic_Guo.md b/docs/_pages/examples/CB_008_anharmonic_Guo.md
index a200d38..cf414fd 100644
--- a/docs/_pages/examples/CB_008_anharmonic_Guo.md
+++ b/docs/_pages/examples/CB_008_anharmonic_Guo.md
@@ -77,4 +77,6 @@ title: ""
     ylim([0,150])
     set(gca,'ydir','reverse')
   end
+
+  saveas(gcf,'./figures/CB_008_anharmonic_Gu0.png')
 ```
diff --git a/docs/_pages/examples/CB_009_anhporo.md b/docs/_pages/examples/CB_009_anhporo.md
index d758afc..0d8c484 100644
--- a/docs/_pages/examples/CB_009_anhporo.md
+++ b/docs/_pages/examples/CB_009_anhporo.md
@@ -4,6 +4,9 @@ title: ""
 ---
 
 # CB_009_anhporo.m
+## output figures
+
+!['CB_009_anhporo'](/vbr/assets/images/CBs/CB_009_anhporo.png){:class="img-responsive"}
 ## contents
 ```matlab
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -42,4 +45,6 @@ title: ""
   semilogx(VBR.in.SV.phi,VBR.out.elastic.anh_poro.Vsu/1e3,'k','linewidth',1.5)
   xlabel('\phi'); ylabel('Vsu(P,T,\phi) [km/s]')
   set(gca,'linewidth',1.5)
+
+  saveas(gcf,'./figures/CB_009_anhporo.png')
 ```
diff --git a/docs/_pages/examples/CB_010_depthProfiles.md b/docs/_pages/examples/CB_010_depthProfiles.md
index 8d231ae..e8e3ee8 100644
--- a/docs/_pages/examples/CB_010_depthProfiles.md
+++ b/docs/_pages/examples/CB_010_depthProfiles.md
@@ -10,7 +10,7 @@ function [VBR,HF] = CB_010_depthProfiles()
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % [VBR,HF] = CB_010_depthProfiles()
   %
-  %   Calculate seismic properties for a half space cooling profile at 80 Myrs
+  %   Calculate seismic properties for a half space cooling profile at 30 Myrs
   %   Assumes constant density for thermal profile (analytical half space
   %   cooling), re-calculates density for the seismic properties.
   %
@@ -58,7 +58,7 @@ function [VBR,HF] = CB_010_depthProfiles()
 
   %% Build figures %%
 
-    figure()
+    figure('PaperPosition',[0,0,8,4],'PaperPositionMode','manual')
     ax1=subplot(1,5,1);
     plot(HF.T_C,HF.z_km)
     hold on
@@ -101,6 +101,7 @@ function [VBR,HF] = CB_010_depthProfiles()
       xlim([0,8])
     end
 
+    saveas(gcf,'./figures/CB_00210_depthProfiles.png')
 end
 
 function HF = HalfspaceModel(age_Myrs)
diff --git a/docs/_pages/examples/CB_011_meltEffects.md b/docs/_pages/examples/CB_011_meltEffects.md
index eae20db..f2a9653 100644
--- a/docs/_pages/examples/CB_011_meltEffects.md
+++ b/docs/_pages/examples/CB_011_meltEffects.md
@@ -4,9 +4,6 @@ title: ""
 ---
 
 # CB_011_meltEffects.m
-## output figures
-
-!['CB_011_meltEffects_case2_fig2'](/vbr/assets/images/CBs/CB_011_meltEffects_case2_fig2.png){:class="img-responsive"}
 ## contents
 ```matlab
 function Results = CB_011_meltEffects(case2run)
diff --git a/docs/_pages/examples/CB_013_G_K_inputs.md b/docs/_pages/examples/CB_013_G_K_inputs.md
new file mode 100644
index 0000000..082df63
--- /dev/null
+++ b/docs/_pages/examples/CB_013_G_K_inputs.md
@@ -0,0 +1,44 @@
+---
+permalink: /examples/CB_013_G_K_inputs/
+title: ""
+---
+
+# CB_013_G_K_inputs.m
+## contents
+```matlab
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% CB_013_G_K_inputs.m
+%
+%  Specify unrelaxed shear and bulk moduli.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% put the VBRc in the path %
+clear
+path_to_top_level_vbr='../../';
+addpath(path_to_top_level_vbr)
+vbr_init
+
+% specify state variables as usual
+VBR = struct();
+VBR.in.SV.T_K = linspace(800, 1000, 4);
+sz_T = size(VBR.in.SV.T_K);
+VBR.in.SV.P_GPa = linspace(2, 3, 4);
+VBR.in.SV.rho = 3300 * ones(sz_T);
+VBR.in.SV.phi = 0.01 * ones(sz_T);
+VBR.in.SV.sig_MPa = 1 * ones(sz_T);
+VBR.in.SV.dg_um = 1e4 * ones(sz_T);
+VBR.in.SV.f = [0.01, 0.1];
+
+% specify methods as usual
+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'};
+
+% also specify the unrelaxed moduli at elevated temperature, pressure.
+% you could, instead, load these from a file!
+VBR.in.elastic.Gu_TP = linspace(50, 60, 4) * 1e9; % shear modulus
+VBR.in.elastic.Ku_TP = VBR.in.elastic.Gu_TP * 1.5; % bulk modulus
+
+% call the VBRc
+VBR = VBR_spine(VBR);
+```
diff --git a/docs/_pages/examples/vbrcore.md b/docs/_pages/examples/vbrcore.md
index 93c5f41..d2e9e27 100644
--- a/docs/_pages/examples/vbrcore.md
+++ b/docs/_pages/examples/vbrcore.md
@@ -17,3 +17,4 @@ Simple "Cookbook" (CB) scripts demonstrating various use cases of the VBR Calcul
 * `CB_010_depthProfiles.m` [link](/vbr/examples/CB_010_depthProfiles/)
 * `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/)
diff --git a/docs/_pages/vbrmethods/el/anharmonic.md b/docs/_pages/vbrmethods/el/anharmonic.md
index fff3ec3..4bc2119 100644
--- a/docs/_pages/vbrmethods/el/anharmonic.md
+++ b/docs/_pages/vbrmethods/el/anharmonic.md
@@ -97,3 +97,20 @@ Gu = Gu_0_ol .* VBR.in.SV.chi + (1-VBR.in.SV.chi) .* Gu_0_crust;
 If `VBR.in.SV.chi` is not set by the user, then  `VBR.in.SV.chi` is initialized to a value of 1 everywhere.
 
 Thus, to produce depth profiles, the state variable arrays should correspond to some depth dependence. For an example, see `Projects/vbr_core_examples/CB_008_anharmonic_Guo.m`.
+
+# Setting unrelaxed moduli at T, P directly
+
+The VBRc relies on a simple linear calculation of the unrelaxed moduli at the
+temperature and pressure of interest. It is possible, however, to load in
+unrelaxed moduli calculated with other programs (like. e.g., Perplex). To do so,
+you can set the following fields:
+
+```matlab
+VBR.in.elastic.Gu_TP = ...
+VBR.in.elastic.Ku_TP = ...
+```
+
+Both `Gu_TP` and `Ku_TP` should be arrays of the same shape as the state variable
+arrays. When these fields are present, the anharmonic calculation will simply read
+from these fields, allowing you to set their values in any way you see fit (e.g.,
+reading from Perplex output or calling your own functions).
diff --git a/docs/assets/images/CBs/CB_001_0D_scalar.png b/docs/assets/images/CBs/CB_001_0D_scalar.png
index 617ae92..ab3b022 100644
Binary files a/docs/assets/images/CBs/CB_001_0D_scalar.png and b/docs/assets/images/CBs/CB_001_0D_scalar.png differ
diff --git a/docs/assets/images/CBs/CB_002_2D_HalfSpaceCooling.png b/docs/assets/images/CBs/CB_002_2D_HalfSpaceCooling.png
new file mode 100644
index 0000000..44e9819
Binary files /dev/null and b/docs/assets/images/CBs/CB_002_2D_HalfSpaceCooling.png differ
diff --git a/docs/assets/images/CBs/CB_003_JF10.png b/docs/assets/images/CBs/CB_003_JF10.png
index f8f9ec3..1fde302 100644
Binary files a/docs/assets/images/CBs/CB_003_JF10.png and b/docs/assets/images/CBs/CB_003_JF10.png differ
diff --git a/docs/assets/images/CBs/CB_004_xfit_premelt.png b/docs/assets/images/CBs/CB_004_xfit_premelt.png
new file mode 100644
index 0000000..5121429
Binary files /dev/null and b/docs/assets/images/CBs/CB_004_xfit_premelt.png differ
diff --git a/docs/assets/images/CBs/CB_005_grainsize_melt.png b/docs/assets/images/CBs/CB_005_grainsize_melt.png
index af0f0ce..616b57f 100644
Binary files a/docs/assets/images/CBs/CB_005_grainsize_melt.png and b/docs/assets/images/CBs/CB_005_grainsize_melt.png differ
diff --git a/docs/assets/images/CBs/CB_006_viscosity.png b/docs/assets/images/CBs/CB_006_viscosity.png
new file mode 100644
index 0000000..6139fe7
Binary files /dev/null and b/docs/assets/images/CBs/CB_006_viscosity.png differ
diff --git a/docs/assets/images/CBs/CB_007_SLB2005.png b/docs/assets/images/CBs/CB_007_SLB2005.png
new file mode 100644
index 0000000..4b13a92
Binary files /dev/null and b/docs/assets/images/CBs/CB_007_SLB2005.png differ
diff --git a/docs/assets/images/CBs/CB_009_anhporo.png b/docs/assets/images/CBs/CB_009_anhporo.png
new file mode 100644
index 0000000..c852c27
Binary files /dev/null and b/docs/assets/images/CBs/CB_009_anhporo.png differ
diff --git a/docs/assets/images/CBs/CB_012_simplecrust.png b/docs/assets/images/CBs/CB_012_simplecrust.png
index fc34eb5..f1c7370 100644
Binary files a/docs/assets/images/CBs/CB_012_simplecrust.png and b/docs/assets/images/CBs/CB_012_simplecrust.png differ