diff --git a/lisainstrument/instrument.py b/lisainstrument/instrument.py
index c13a260d17ca436c0940672b1bd29c4adaa3e0b2..2928d0853c1633344d610a446e7f5f5ea564714c 100644
--- a/lisainstrument/instrument.py
+++ b/lisainstrument/instrument.py
@@ -34,7 +34,7 @@ class Instrument:
     SCS = ForEachSC.indices()
     MOSAS = ForEachMOSA.indices()
 
-    def __init__(self, size=2592000, dt=1/4, t0=0,
+    def __init__(self, size=2592000, dt=1/4, t0='orbits',
                  # Inter-spacecraft propagation
                  orbits=None, gws=None, glitches=None, interpolation=('lagrange', 31),
                  # Lasers
@@ -54,7 +54,7 @@ class Instrument:
         Args:
             size: number of samples to generate
             dt: sampling period [s]
-            t0: initial time [s]
+            t0: initial time [s], or 'orbits' to match that of the orbits
             orbits: path to orbit file, or dictionary of scalars or time series for PPRs,
                 or None for default set of PPRs corresponding to static arms fit from Keplerian orbits
                 (from LISA Orbits v1.0) around t = 0
@@ -89,7 +89,7 @@ class Instrument:
                 ('kaiser', attenuation [dB], f1 [Hz], f2 [Hz]) with f1 < f2 the frequencies defining
                 the transition band
         """
-        # pylint: disable=too-many-arguments,too-many-statements
+        # pylint: disable=too-many-arguments,too-many-statements,too-many-locals,too-many-branches
         logger.info("Initializing instrumental simulation")
         self.git_url = 'https://gitlab.in2p3.fr/lisa-simulation/instrument'
         self.version = meta.__version__
@@ -98,7 +98,15 @@ class Instrument:
         # Measurement sampling
         self.size = int(size)
         self.dt = float(dt)
-        self.t0 = float(t0)
+        if t0 == 'orbits':
+            if isinstance(orbits, str):
+                logger.debug("Reading initial time from orbit file '%s'", orbits)
+                with h5py.File(orbits, 'r') as orbitf:
+                    self.t0 = float(orbitf.attrs['tau0'])
+            else:
+                self.t0 = 0.0
+        else:
+            self.t0 = float(t0)
         self.fs = 1 / self.dt
         logger.info("Computing measurement time vector (size=%s, dt=%s)", self.size, self.dt)
         self.t = self.t0 + numpy.arange(self.size, dtype=numpy.float64) * self.dt