diff --git a/DESCRIPTION b/DESCRIPTION
index 3b12abb..ccb7d25 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: VAST
 Type: Package
 Title: Vector-Autoregressive Spatio-Temporal (VAST) Model
-Version: 3.10.0
+Version: 3.10.1
 Date: 2022-11-18
 Authors@R: 
     c(person(given = "James",
@@ -30,7 +30,7 @@ Imports:
     devtools
 Depends:
     TMB (>= 1.8.0),
-    FishStatsUtils,
+    FishStatsUtils (>= 2.12.1),
     R (>= 3.5.0)
 Suggests:
     testthat,
diff --git a/R/apply_epsilon.R b/R/apply_epsilon.R
index d401acb..0382be3 100644
--- a/R/apply_epsilon.R
+++ b/R/apply_epsilon.R
@@ -36,6 +36,7 @@ function( fit,
   Data = Obj$env$data
   Map = Obj$env$map
   Random = fit$tmb_list$Random
+  Version_framework = paste0( fit$input_args$model_args_input$Version, "_", fit$input_args$model_args_input$framework )
 
   # Extract and modify parameters
   New_params = fit$ParHat
@@ -53,7 +54,8 @@ function( fit,
                     map = Map,
                     random = Random,
                     intern = TRUE,
-                    inner.control = inner.control )
+                    inner.control = inner.control,
+                    DLL = Version_framework )
   obj$env$beSilent()
   gradient = obj$gr(fixed)
 
diff --git a/R/make_data.R b/R/make_data.R
index 3253133..d2f09ea 100644
--- a/R/make_data.R
+++ b/R/make_data.R
@@ -859,7 +859,7 @@ function( b_i,
   # Mirroring
   if( RhoConfig[4]==6 ){
     if( FieldConfig_input[2,1]!=FieldConfig_input[2,2] ){
-      stop("To fix 'Epsilon_rho2_f` equal to 'Epsilon_rho2_f`, you must specify the same rank using `FieldConfig_input[2,1]` equal to `FieldConfig_input[2,2]`")
+      stop("To fix 'Epsilon_rho2_f` equal to 'Epsilon_rho1_f`, you must specify the same rank using `FieldConfig_input[2,1]` equal to `FieldConfig_input[2,2]`")
     }
   }
 
diff --git a/R/project_model.R b/R/project_model.R
index feaa5df..80a5d85 100644
--- a/R/project_model.R
+++ b/R/project_model.R
@@ -86,6 +86,14 @@ function( x,
     }
   }
 
+  # Errors
+  if( any(x$data_list$RhoConfig %in% c(0,3)) ){
+    #stop("`project_model` is currently designed to work only with temporally varying epsilon and beta terms")
+  }
+  if( any(x$data_list$RhoConfig[c("Beta1","Beta2")] %in% c(0)) ){
+    stop("`project_model` is currently designed to work only with temporally varying or constant beta terms")
+  }
+
   ##############
   # Step 1: Generate uncertainty in historical period
   ##############
@@ -125,15 +133,16 @@ function( x,
   # Step 2: Generate uncertainty in historical period
   ##############
 
-  t_i = c( x$data_frame$t_i, max(x$data_frame$t_i)+rep(1:n_proj,each=2) )
+  t_i = c( x$data_frame$t_i, max(x$data_frame$t_i)+rep(seq_len(n_proj),each=2) )
   b_i = c( x$data_list$b_i, as_units(rep(c(0,mean(x$data_frame$b_i)),n_proj), units(x$data_list$b_i)) )
   v_i = c( x$data_frame$v_i, rep(0,2*n_proj) )
   Lon_i = c( x$data_frame$Lon_i, rep(mean(x$data_frame$Lon_i),2*n_proj) )
   Lat_i = c( x$data_frame$Lat_i, rep(mean(x$data_frame$Lat_i),2*n_proj) )
   a_i = c( x$data_list$a_i, as_units(rep(mean(x$data_frame$a_i),2*n_proj), units(x$data_list$a_i)) )
   PredTF_i = c( x$data_list$PredTF_i, rep(1,2*n_proj) )
-  c_iz = rbind( x$data_list$c_iz, x$data_list$c_iz[rep(1:n_proj,each=2),,drop=FALSE] )
-  new_catchability_data = rbind( x$catchability_data, x$catchability_data[rep(1:n_proj,each=2),,drop=FALSE] )
+  c_iz = rbind( x$data_list$c_iz, x$data_list$c_iz[rep(seq_len(n_proj),each=2),,drop=FALSE] )
+  new_catchability_data = rbind( x$catchability_data, x$catchability_data[rep(seq_len(n_proj),each=2),,drop=FALSE] )
+  proj_t = x$data_list$n_t + seq_len( n_proj )
 
   ##############
   # Step 3: Build object with padded bounds
@@ -190,6 +199,18 @@ function( x,
       }
     }
 
+    # Deal with beta1/beta2 = 3
+    if( x$data_list$RhoConfig["Beta1"]==3 ){
+      tmp = ParList1$beta1_ft
+      tmp[,proj_t] = NA
+      ParList1$beta1_ft = ifelse( is.na(tmp), rowMeans(tmp,na.rm=TRUE)%o%rep(1,ncol(tmp)), ParList1$beta1_ft )
+    }
+    if( x$data_list$RhoConfig["Beta2"]==3 ){
+      tmp = ParList1$beta2_ft
+      tmp[,proj_t] = NA
+      ParList1$beta2_ft = ifelse( is.na(tmp), rowMeans(tmp,na.rm=TRUE)%o%rep(1,ncol(tmp)), ParList1$beta2_ft )
+    }
+
     ##############
     # Step 5: Re-build model
     ##############
@@ -217,6 +238,25 @@ function( x,
       Parameters = ParList1,
       working_dir = working_dir )
 
+    # Debugging comparisons
+    if( FALSE ){
+      # Samples
+      range(u_zr[,1] - Obj$env$last.par.best)
+      # Compare fixed effects
+      v0 - x$tmb_list$Obj$env$last.par.best[-x$tmb_list$Obj$env$random]
+      v2 = x2$tmb_list$Obj$env$last.par.best[-x2$tmb_list$Obj$env$random]
+      range(v0 - v2)
+      # Compare parameters
+      v0 = unlist(x$tmb_list$Parameters[-match(x$tmb_list$Random,names(ParList))])
+      #v1 = unlist(ParList1[-match(x1$tmb_list$Random,names(ParList))])
+      v2 = unlist(x2$tmb_list$Parameters[-match(x2$tmb_list$Random,names(ParList))])
+      range(v1-v2)
+      # compare last.par.best
+      v0 = x$tmb_list$Obj$env$last.par.best[-x$tmb_list$Obj$env$random]
+      v2 = x2$tmb_list$Obj$env$last.par.best[-x2$tmb_list$Obj$env$random]
+      range( v0 - v2 )
+    }
+
     ##############
     # Step 5: Simulate random effects
     ##############
diff --git a/inst/executables/VAST_v14_0_1.cpp b/inst/executables/VAST_v14_0_1.cpp
index a4f18b4..184ec04 100644
--- a/inst/executables/VAST_v14_0_1.cpp
+++ b/inst/executables/VAST_v14_0_1.cpp
@@ -210,9 +210,10 @@ Type rgengamma( Type mean,
   // See: C:\Users\James.Thorson\Desktop\Work files\AFSC\2021-10 -- Generalized gamma-lognormal\Explore gengamma.R
   Type k = pow( lambda, -2 );
   Type Shape = pow( sigma, -1 ) * lambda;
-  Type Scale = mean / exp(lgamma( (k*Shape+1)/Shape )) * exp(lgamma( k ));
+  //Type Scale = mean / exp(lgamma( (k*Shape+1)/Shape )) * exp(lgamma( k ));
+  Type log_Scale = log(mean) - lgamma( (k*Shape+1)/Shape ) + lgamma( k );
   Type w = log(rgamma(k, Type(1.0)));
-  Type y = w/Shape + log(Scale);
+  Type y = w/Shape + log_Scale;
   return exp(y);
 }
 
@@ -1648,8 +1649,10 @@ Type objective_function<Type>::operator() ()
       jnll_beta1 -= dnorm( beta1_tf(t,f), beta1_mean_tf(t,f), Type(1.0), true );
       // Simulate new values when using obj.simulate()
       if( (Options(14) == 1) | (simulate_t(t) == 1) ){
-        SIMULATE{
-          beta1_tf(t,f) = rnorm( beta1_mean_tf(t,f), Type(1.0) );
+        if( (RhoConfig(0)==1) | (RhoConfig(0)==2) | (RhoConfig(0)==4) | (RhoConfig(0)==5) ){
+          SIMULATE{
+            beta1_tf(t,f) = rnorm( beta1_mean_tf(t,f), Type(1.0) );
+          }
         }
       }
     }}
@@ -1691,8 +1694,10 @@ Type objective_function<Type>::operator() ()
       jnll_beta2 -= dnorm( beta2_tf(t,f), beta2_mean_tf(t,f), Type(1.0), true );
       // Simulate new values when using obj.simulate()
       if( (Options(14) == 1) | (simulate_t(t) == 1) ){
-        SIMULATE{
-          beta2_tf(t,f) = rnorm( beta2_mean_tf(t,f), Type(1.0) );
+        if( (RhoConfig(1)==1) | (RhoConfig(1)==2) | (RhoConfig(1)==4) | (RhoConfig(1)==5) | (RhoConfig(1)==6) ){
+          SIMULATE{
+            beta2_tf(t,f) = rnorm( beta2_mean_tf(t,f), Type(1.0) );
+          }
         }
       }
     }}
@@ -2032,7 +2037,7 @@ Type objective_function<Type>::operator() ()
           }
           // Generalized-gamma;  mean, sigma, lambda parameterization
           if(ObsModel_ez(e_i(i),0)==9){
-            LogProb2_i(i) = dgengamma(b_i(i), R2_i(i), SigmaM(e_i(i),0), logSigmaM(e_i(i),1), true);
+            LogProb2_i(i) = dgengamma(b_i(i), R2_i(i), SigmaM(e_i(i),0), SigmaM(e_i(i),1), true);
             deviance2_i(i) = NAN;
             // Simulate new values when using obj.simulate()
             // Could be updated, available as rgengamma.orig
diff --git a/manual/NEWS.docx b/manual/NEWS.docx
index 01f3096..41ab645 100644
Binary files a/manual/NEWS.docx and b/manual/NEWS.docx differ