diff --git a/Exec/science/subchandra/_prob_params b/Exec/science/subchandra/_prob_params
index 627ea5c7ea..e546592e9f 100644
--- a/Exec/science/subchandra/_prob_params
+++ b/Exec/science/subchandra/_prob_params
@@ -18,3 +18,13 @@ mu_p                  real           -5.0_rt         y
 mu_n                  real           -11.0_rt        y
 
 ihe4                  integer         -1
+
+# second perturbation angle -- if it is negative, then it is disabled
+# this is only relevant in 3D
+second_pert_angle         real           -1              y
+
+second_R_pert             real           1.e7_rt         y
+
+second_pert_temp_factor   real           10.0_rt         y
+
+second_pert_rad_factor    real           2.0_rt          y
diff --git a/Exec/science/subchandra/problem_initialize_state_data.H b/Exec/science/subchandra/problem_initialize_state_data.H
index 249d151931..95a1a7178d 100644
--- a/Exec/science/subchandra/problem_initialize_state_data.H
+++ b/Exec/science/subchandra/problem_initialize_state_data.H
@@ -70,7 +70,7 @@ void problem_initialize_state_data (int i, int j, int k,
     state(i,j,k,UMY) = 0.0_rt;
     state(i,j,k,UMZ) = 0.0_rt;
 
-    // add a perturbation
+    // add a perturbation at the north pole
 
     Real t0 = state(i,j,k,UTEMP);
 
@@ -96,6 +96,36 @@ void problem_initialize_state_data (int i, int j, int k,
                                (0.150e0_rt * (1.0_rt + std::tanh(2.0_rt - r1))));
 
 
+#if AMREX_SPACEDIM == 3
+    // optional second perturbation
+
+    if (problem::second_pert_angle > 0) {
+
+        Real angle_rad = problem::second_pert_angle * M_PI / 180.0_rt;
+
+        Real second_pert_center = problem::second_R_pert + problem::R_He_base;
+
+        // these are measured with respect to the center
+
+        Real x_pert = second_pert_center * std::sin(angle_rad);
+        Real y_pert = 0.0_rt;
+        Real z_pert = second_pert_center * std::cos(angle_rad);
+
+        Real dx_pert = x - x_pert;
+        Real dy_pert = y - y_pert;
+        Real dz_pert = z - z_pert;
+
+        Real r2 = std::sqrt(dx_pert * dx_pert + dy_pert * dy_pert + dz_pert * dz_pert) /
+            (2.5e6_rt * problem::second_pert_rad_factor);
+
+        // we are assuming here that the 2 perturbations don't overlap
+
+        state(i,j,k,UTEMP) = t0 * (1.0_rt + X_he * problem::second_pert_temp_factor *
+                                   (0.150e0_rt * (1.0_rt + std::tanh(2.0_rt - r2))));
+
+    }
+#endif
+
     burn_state.rho = state(i,j,k,URHO);
     burn_state.T = state(i,j,k,UTEMP);
     // we don't need to refill xn, since it still holds unchanged from above