diff --git a/flopy/modflow/mfhob.py b/flopy/modflow/mfhob.py
index 5d0a9a89a8..19f1cae593 100755
--- a/flopy/modflow/mfhob.py
+++ b/flopy/modflow/mfhob.py
@@ -371,7 +371,8 @@ def load(cls, f, model, ext_unit_dict=None, check=True):
         nobs = 0
 
         # set to False for 1st call to ensure that totim cache is updated
-        use_cached_totim = False
+        tmax = model.dis.get_final_totim()
+        use_cached_totim = True
 
         while True:
             # read dataset 3
@@ -431,7 +432,6 @@ def load(cls, f, model, ext_unit_dict=None, check=True):
                 names = [obsnam]
                 tsd = [totim, hob]
                 nobs += 1
-                use_cached_totim = True
             else:
                 names = []
                 tsd = []
@@ -452,7 +452,6 @@ def load(cls, f, model, ext_unit_dict=None, check=True):
                     hob = float(t[3])
                     tsd.append([totim, hob])
                     nobs += 1
-                    use_cached_totim = True
 
             obs_data.append(
                 HeadObservation(
@@ -466,6 +465,7 @@ def load(cls, f, model, ext_unit_dict=None, check=True):
                     obsname=obsnam,
                     mlay=mlay,
                     itt=itt,
+                    tmax=tmax,
                     time_series_data=tsd,
                     names=names,
                 )
@@ -543,6 +543,9 @@ class HeadObservation:
         observations. itt = 1 specified for heads and itt = 2 specified
         if initial value is head and subsequent changes in head. Only
         specified if irefsp is < 0. Default is 1.
+    tmax : float
+        Maximum simulation time calculated using get_final_totim function of 
+        ModflowDis. Added to avoid repetitive calls. 
     mlay : dictionary of length (abs(irefsp))
         Key represents zero-based layer numbers for multilayer observations and
         value represents the fractional value for each layer of multilayer
@@ -588,6 +591,7 @@ def __init__(
         roff=0.0,
         coff=0.0,
         itt=1,
+        tmax=None,
         mlay=None,
         time_series_data=None,
         names=None,
@@ -645,7 +649,8 @@ def __init__(
             time_series_data = np.reshape(time_series_data, (1, 2))
 
         # find indices of time series data that are valid
-        tmax = model.dis.get_final_totim()
+        if tmax is None:
+            tmax = model.dis.get_final_totim()
         keep_idx = time_series_data[:, 0] <= tmax
         time_series_data = time_series_data[keep_idx, :]
 
@@ -677,22 +682,18 @@ def __init__(
                     "{} names are required.".format(len(names), self.nobs)
                 )
 
-        # set use_cached_totim to False first to ensure totim is updated
-        use_cached_totim = False
-
         # create time_series_data
         self.time_series_data = self._get_empty(ncells=shape[0])
         for idx in range(self.nobs):
             t = time_series_data[idx, 0]
             kstp, kper, toffset = model.dis.get_kstp_kper_toffset(
-                t, use_cached_totim
+                t, use_cached_totim=True
             )
             self.time_series_data[idx]["totim"] = t
             self.time_series_data[idx]["irefsp"] = kper
             self.time_series_data[idx]["toffset"] = toffset / tomulth
             self.time_series_data[idx]["hobs"] = time_series_data[idx, 1]
             self.time_series_data[idx]["obsname"] = names[idx]
-            use_cached_totim = True
 
         if self.nobs > 1:
             self.irefsp = -self.nobs