diff --git a/R/VAST_do.R b/R/VAST_do.R
index da47721..5b60167 100644
--- a/R/VAST_do.R
+++ b/R/VAST_do.R
@@ -74,6 +74,8 @@ VAST_do <- function(Database, settings, conditiondir, compiledir,
     fine_scale = settings[["fine_scale"]],
     strata.limits = settings[["strata"]],
     #zone = NA, #default
+    # automatically detecting zone in `FishStatsUtils::make_extrapolation_info()`
+    # DO NOT allow users to control just use default
     FieldConfig = settings[["FieldConfig"]],
     RhoConfig = settings[["RhoConfig"]],
     OverdispersionConfig = settings[["overdispersion"]],
@@ -81,7 +83,7 @@ VAST_do <- function(Database, settings, conditiondir, compiledir,
     bias.correct = TRUE,
     #calculate derived quantities of interest, no harm
     #Options = ,#default is more than original VASTWestCoast
-    #use_anisotropy = TRUE, #default
+    use_anisotropy = settings[["use_anisotropy"]],
     Version = settings[["version"]],
     #treat_nonencounter_as_zero = TRUE, #default
     #VamConfig = c(Method = 0, Rank = 0, Timing = 0), #default
@@ -130,7 +132,8 @@ VAST_do <- function(Database, settings, conditiondir, compiledir,
     model_args = list(CompileDir = compiledir),
     silent = TRUE,
     run_model = TRUE), error = function(e) e)
-  if ("simpleError" %in% class(out)) {
+  if ("simpleError" %in% class(out) |
+      is.null(out$parameter_estimates$SD)) {
     save(out, file = file.path(conditiondir, "Error.RData"))
     save(list = ls(all.names = TRUE), file = file.path(conditiondir, "Save.RData"))
     return(out)
@@ -142,6 +145,7 @@ VAST_do <- function(Database, settings, conditiondir, compiledir,
     error = function(e) e)
 
   index <- suppressWarnings(FishStatsUtils::plot_biomass_index(
+    fit = out,
     DirName = file.path(conditiondir, .Platform$file.sep),
     TmbData = out$data_list, Sdreport = out$parameter_estimates$SD,
     use_biascorr = TRUE,
@@ -154,8 +158,16 @@ VAST_do <- function(Database, settings, conditiondir, compiledir,
     parameter_estimates = out$parameter_estimates,
     savedir = conditiondir)
   fileindex <- file.path(conditiondir, "Table_for_SS3.csv")
-  indexdata <- utils::read.csv(fileindex)
-  indexdata[,"Year"] <- out$year_labels[indexdata[,"Year"]]
+  indexdata <- data.frame(
+    Year = out$year_labels[index[["Table"]][, "Time"]],
+    # Change units to mt since converting from kg to mt below
+    Unit = "mt",
+    Fleet = index[["Table"]][, "Stratum"],
+    # Go from kg to mt to keep backwards compatibility with VASTWestCoast
+    Estimate_metric_tons = index[["Table"]][, "Estimate"] / 1000,
+    SD_log = index[["Table"]][, "Std. Error for ln(Estimate)"],
+    SD_mt = index[["Table"]][, "Std. Error for Estimate"] / 1000
+  )
   utils::write.csv(x = indexdata, file = fileindex, row.names = FALSE)
   plot_ss(file.in = fileindex,
     lab.survey = survey,
diff --git a/R/VAST_spp.R b/R/VAST_spp.R
index 21a6c45..039fd00 100644
--- a/R/VAST_spp.R
+++ b/R/VAST_spp.R
@@ -85,7 +85,7 @@
 #' }
 #'
 VAST_spp <- function(dir, species,
-  dist = c("lognormal", "gamma")) {
+  dist = c("lognormal", "gamma", "tweedie")) {
 
   #### Get species, survey, and strata info
   info <- nwfscSurvey::GetSpp.fn(species)
@@ -107,7 +107,8 @@ VAST_spp <- function(dir, species,
   for (obs in dist) {
     obs_model <- switch(obs,
       lognormal = c(1, 0),
-      gamma = c(2, 0)
+      gamma = c(2, 0),
+      tweedie = c(10, 2)
     )
     modeldir <- file.path(sppdir, paste(survey, obs, sep = "_"))
     settings <- list(
diff --git a/R/get_settings.R b/R/get_settings.R
index e7e1203..6e67c5a 100644
--- a/R/get_settings.R
+++ b/R/get_settings.R
@@ -32,7 +32,9 @@ get_settings <- function(settings = NULL, verbose = FALSE) {
     "field" = NULL,
     "rho" = NULL,
     "fine_scale" = TRUE,
-    "overdispersion" = NULL)
+    "overdispersion" = NULL,
+    "use_anisotropy" = TRUE
+    )
   need <- !names(Settings_all) %in% names(settings)
   if (verbose) {
     message("Adding the following objects to settings:\n",
@@ -67,6 +69,11 @@ get_settings <- function(settings = NULL, verbose = FALSE) {
     default = c(Beta1 = 0, Beta2 = 0, Epsilon1 = 0, Epsilon2 = 0))
   Settings_all$FieldConfig <- get_settings_single(Settings_all$field,
     default = c(Omega1 = 1, Epsilon1 = 1, Omega2 = 1, Epsilon2 = 1))
+  if (Settings_all[["ObsModelcondition"]][1] == 10) {
+    # Following items are specific to the Tweedie distribution #43
+    Settings_all[["FieldConfig"]][c("Omega1", "Epsilon1")] <- 0
+    Settings_all[["RhoConfig"]][c("Beta1")] <- 3
+  }
 
   # Overdispersion
   if (is.null(Settings_all[["overdispersion"]])) {
diff --git a/R/plot_qq.R b/R/plot_qq.R
index 49803a2..33c922f 100644
--- a/R/plot_qq.R
+++ b/R/plot_qq.R
@@ -24,24 +24,28 @@ plot_qq <- function(outdir, model) {
   n_e <- model[["data_list"]][["n_c"]]
   e_i <- model[["data_list"]][["n_e"]]
   sigmaM <- (rep(1, n_e) %o% model[["Report"]][["SigmaM"]])[1, 1, 1]
-  Which <- which(model[["data_list"]][["b_i"]] > 0)
+  Which <- which(as.vector(model[["data_list"]][["b_i"]]) > 0)
 
   if (model[["data_list"]][["ObsModel_ez"]][1] == 1) {
     Qvals <- na.omit(mapply(plnorm,
-          q = model[["data_list"]][["b_i"]][Which],
-          meanlog = log(model[["data_list"]][["a_i"]][Which] * exp(model[["Report"]][["P2_iz"]][Which])) -
+          q = as.vector(model[["data_list"]][["b_i"]])[Which],
+          meanlog = log(as.vector(model[["data_list"]][["a_i"]])[Which] * exp(model[["Report"]][["P2_iz"]][Which])) -
                     (sigmaM^2) / 2,
           MoreArgs = list(sdlog = sigmaM)))
   }
   if (model[["data_list"]][["ObsModel_ez"]][1] == 2) {
     Qvals <- na.omit(mapply(pgamma,
-          q = model[["data_list"]][["b_i"]][Which],
+          q = as.vector(model[["data_list"]][["b_i"]])[Which],
           scale = (sigmaM^2) *
-            (model[["data_list"]][["a_i"]][Which] *
+            (as.vector(model[["data_list"]][["a_i"]])[Which] *
             exp(model[["Report"]][["P2_iz"]][Which])
             ),
           MoreArgs = list(shape = 1 / (sigmaM^2))))
   }
+  if (model[["data_list"]][["ObsModel_ez"]][1] == 10) {
+    message("QQ plot is not yet implemented for the Tweedie distribution. See #45.")
+    return(invisible())
+  }
 
   grDevices::png(file.path(outdir, "VASTWestCoast_QQ.png"), height = 18, width = 18, units = "cm",
     res = 300)
diff --git a/R/plot_ss.R b/R/plot_ss.R
index b940f5f..6a8c8d0 100644
--- a/R/plot_ss.R
+++ b/R/plot_ss.R
@@ -49,19 +49,19 @@ plot_ss <- function(file.in = "Table_for_SS3.csv", savefile,
   savefile <- normalizePath(savefile, mustWork = FALSE)
   dir.create(dirname(savefile), recursive = TRUE, showWarnings = FALSE)
 
-  colors <- c("black",  "blue", "darkorchid1", "red")
+  colors <- c("black",  "red", "darkorchid1", "blue")
   pch.vec <- c(21, 22, 23, 24)
   cex.vec <- c(1.6, 1.4, 1.4, 1.4)
   cex.lab  <- 2.1
   cex.axis <- 1.6
-  colors <- c("black",  "blue", "darkorchid1", "red")
+
   if (length(strat) > 4) {
     stop("plot_ss is not designed to work with more than 4 strata.",
       "\nPlease contact the maintainer(s) to request additional",
       "\nfunctionality here or limit the number of strata in the csv.")
   }
   if ("north" %in% tolower(strat) && length(strat) == 3) {
-    colors <- c("black", "blue", "red")
+    colors <- c("black", "red", "blue")
   }
 
   grDevices::png(filename = savefile,
diff --git a/R/summary_nwfsc.R b/R/summary_nwfsc.R
index d2f8937..9f81cb4 100644
--- a/R/summary_nwfsc.R
+++ b/R/summary_nwfsc.R
@@ -42,7 +42,9 @@ summary_nwfsc <- function(obj, parameter_estimates, savedir = NULL) {
     "Poisson-link", 
     switch(as.character(obj$env$data$ObsModel[1]),
       "1" = "Lognormal",
-      "2" = "Gamma")))
+      "2" = "Gamma",
+      "10" = "Tweedie"))
+  )
   TableA[6, ] <- c("Spatial effect for encounter probability", 
     ifelse(obj$env$data$FieldConfig[1, 2] < 0, "No", "Yes"))
   TableA[7, ] <- c("Spatio-temporal effect for encounter probability", 
diff --git a/man/VAST_spp.Rd b/man/VAST_spp.Rd
index 8e187a5..cce5590 100644
--- a/man/VAST_spp.Rd
+++ b/man/VAST_spp.Rd
@@ -4,7 +4,7 @@
 \alias{VAST_spp}
 \title{Run \pkg{VAST} models for a species sampled by west coast surveys}
 \usage{
-VAST_spp(dir, species, dist = c("lognormal", "gamma"))
+VAST_spp(dir, species, dist = c("lognormal", "gamma", "tweedie"))
 }
 \arguments{
 \item{dir}{An existing directory where you want the new folder that is created