diff --git a/CHANGELOG.md b/CHANGELOG.md
index 66487124..0ec1ff87 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -57,3 +57,5 @@ A collection of bugfixes and feature improvement (big thanks to @evcano for majo
 - output_optim.txt was not writing the correct misfit values for each 
   iteration, fixed
 
+## v2.3.1
+- Hotfix for updated SPECFEM2D parameter file #197
diff --git a/docs/conf.py b/docs/conf.py
index a04be4d8..65cf1369 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -19,9 +19,9 @@
 # -- Project information -----------------------------------------------------
 
 project = 'SeisFlows'
-copyright = '2023'
+copyright = '2024'
 author = 'adjTomo Dev Team'
-version = '2.3.0'
+version = '2.3.1'
 
 
 # -- General configuration ---------------------------------------------------
diff --git a/pyproject.toml b/pyproject.toml
index 9bec584c..8f6511c6 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "seisflows"
-version = "2.3.0"
+version = "2.3.1"
 description = "An automated workflow tool for full waveform inversion"
 readme = "README.md"
 requires-python = ">=3.7"
diff --git a/seisflows/seisflows.py b/seisflows/seisflows.py
index 0a4ed016..e6e08799 100755
--- a/seisflows/seisflows.py
+++ b/seisflows/seisflows.py
@@ -841,7 +841,7 @@ def sempar(self, parameter, value=None, skip_print=False,
         # formatted the same as the rest of the file
         if parameter == "VELOCITY_MODEL":
             key = parameter
-            items = getpar_vel_model(file=par_file)
+            items = getpar_vel_model(file=par_file, strip=True)
             cur_val = ""
         else:
             try:
@@ -861,7 +861,7 @@ def sempar(self, parameter, value=None, skip_print=False,
                 setpar_vel_model(file=par_file, model=value.split("+"))
                 if not skip_print:
                     items.append("->")
-                    items += getpar_vel_model(file=par_file)
+                    items += getpar_vel_model(file=par_file, strip=True)
                     print(msg.cli(f"{key}:", items=items))
             else:
                 setpar(key=parameter, val=value, file=par_file, delim="=")
diff --git a/seisflows/tools/specfem.py b/seisflows/tools/specfem.py
index ae047fb2..814cbda2 100644
--- a/seisflows/tools/specfem.py
+++ b/seisflows/tools/specfem.py
@@ -155,38 +155,69 @@ def setpar(key, val, file, delim="=", match_partial=False):
         f.writelines(lines)
 
 
-def getpar_vel_model(file):
+def _getidx_vel_model(lines):
     """
-    SPECFEM2D doesn't follow standard formatting when defining its internal
-    velocity models so we need a special function to address this specifically.
-    Velocity models are ASSUMED to be formatted in the following way in the
-    SPECFEM2D Par_file (with any number of comment lines in between)
+    Get the line indices of a velocity model, which can be used to retrieve 
+    or replace the model values in a SPECFEM2D paramter file. Used by 
+    `getpar_vel_model` and `setpar_vel_model`
+
+    :type lines: list
+    :param lines: list of strings read from the par_file
+    :rtype idxs: list
+    :param idxs: list of integer indices of the velocity model lines
+    """
+    idxs = []
+    for l, line in enumerate(lines):
+        # Skip over all other parameters, comments, newlines etc.
+        if "=" in line:
+            continue
+        elif line.startswith(" "):
+            continue
+        elif line.startswith("#"):
+            continue
+        elif line.startswith("\n"):
+            continue
+
+        # Matches formatting expected by velocity model and starts with integer
+        # Should be enough to avoid matching other strings
+        lines = line.strip().split()
+        if len(lines) == 15 and lines[0].isdigit():
+            idxs.append(l)
+
+    return idxs
+
+
+def getpar_vel_model(file, strip=False):
+    """
+    SPECFEM2D doesn't follow standard key = val formatting when defining its 
+    internal velocity models so we need a special function to address this 
+    specifically.
+    
+    Velocity models are ASSUMED to be formatted in the following way
 
-    nbmodels                        = 4
     1 1 2700.d0 3000.d0 1732.051d0 0 0 9999 9999 0 0 0 0 0 0
-    2 1 2500.d0 2700.d0 0 0 0 9999 9999 0 0 0 0 0 0
-    3 1 2200.d0 2500.d0 1443.375d0 0 0 9999 9999 0 0 0 0 0 0
-    4 1 2200.d0 2200.d0 1343.375d0 0 0 9999 9999 0 0 0 0 0 0
-    TOMOGRAPHY_FILE                 = ./DATA/tomo_file.xyz
 
+    That is, 15 entries separated by spaces. We use that to find all relevant 
+    lines of the model.
+    
     :type file: str
     :param file: The SPECFEM Par_file to match against
+    :type strip: bool
+    :param strip: strip newline '\n' from each of the model lines
     :rtype: list of str
     :return: list of all the layers of the velocity model as strings
     """
-    _, _, i_start = getpar("nbmodels", file)
-    _, _, i_end = getpar("tomography_file", file)
-
     # i_start + 1 to skip over the 'nbmodels' parameter
-    lines = open(file, "r").readlines()[i_start + 1:i_end]
+    lines = open(file, "r").readlines()
+    idxs = _getidx_vel_model(lines)
     vel_model = []
-    for line in lines:
-        # Skip comments, empty lines, newlines
-        for not_allowed in [" ", "#", "\n"]:
-            if line.startswith(not_allowed):
-                break
+    for idx in idxs:
+        if strip:
+            line = lines[idx].strip()
         else:
-            vel_model.append(line.strip())
+            line = lines[idx]
+        vel_model.append(line)
+
     return vel_model
 
 
@@ -207,34 +238,19 @@ def setpar_vel_model(file, model):
         model = ["1 1 2700.d0 3000.d0 1732.051d0 0 0 9999 9999 0 0 0 0 0 0",
                  "2 1 2500.d0 2700.d0 0 0 0 9999 9999 0 0 0 0 0 0"]
     """
-    _, nbmodels, i_start = getpar("nbmodels", file)
-    i_start += 1  # increase by one to start AFTER nbmodels line
-    _, _, i_end = getpar("tomography_file", file)
-
     lines = open(file, "r").readlines()
-    model_lines = []
-    # i_start + 1 to skip over the 'nbmodels' parameter
-    for i, line in enumerate(lines[i_start:i_end]):
-        # Skip comments, empty lines, newlines
-        for not_allowed in [" ", "#", "\n"]:
-            if line.startswith(not_allowed):
-                break
-        else:
-            # We will use these indices to delete the original model
-            model_lines.append(i)
-
-    # Make sure that our indices are relative to the list and not enumeration
-    model_lines = [_ + i_start for _ in model_lines]
+    model_lines = _getidx_vel_model(open(file, "r").readlines())
 
-    # one-liner to get rid of the original model
+    # one-liner to get rid of the original model from the file
     lines = [i for j, i in enumerate(lines) if j not in model_lines]
+    model_idx_start = model_lines[0]
 
     # Throw a new line onto the last line of the model to get proper formatting
     model[-1] = model[-1] + "\n"
 
     # Drop in the new model one line at a time
     for i, val in enumerate(model):
-        lines.insert(i + model_lines[0], f"{val.strip()}\n")
+        lines.insert(model_idx_start + i, f"{val.strip()}\n")
 
     # Overwrite with new lines containing updated velocity model
     with open(file, "w") as f: