diff --git a/cxx/include/mspass/seismic/CoreSeismogram.h b/cxx/include/mspass/seismic/CoreSeismogram.h
index 247414d67a..685df5359e 100644
--- a/cxx/include/mspass/seismic/CoreSeismogram.h
+++ b/cxx/include/mspass/seismic/CoreSeismogram.h
@@ -223,8 +223,8 @@ class CoreSeismogram : public mspass::seismic::BasicTimeSeries,
   the ends of the rhs are inside the range of the lhs.
 
   \param d is other signal to add to this.
-  \exception MsPASSError can be thrown if lhs and rhs do not have matching
-  time standards.
+  \exception MsPASSError if lhs and rhs have incompatible time standards,
+  sample intervals, or start-time grids.
   **/
   CoreSeismogram &operator+=(const CoreSeismogram &d);
   /*! Addition operator.
@@ -250,8 +250,8 @@ class CoreSeismogram : public mspass::seismic::BasicTimeSeries,
   the ends of the rhs are inside the range of the lhs.
 
   \param d is other signal to subract from this.
-  \exception MsPASSError can be thrown if lhs and rhs do not have matching
-  time standards.
+  \exception MsPASSError if lhs and rhs have incompatible time standards,
+  sample intervals, or start-time grids.
   **/
   CoreSeismogram &operator-=(const CoreSeismogram &d);
   /*! Subtraction operator.
diff --git a/cxx/include/mspass/seismic/CoreTimeSeries.h b/cxx/include/mspass/seismic/CoreTimeSeries.h
index b96b2dd25e..5a3a6d792a 100644
--- a/cxx/include/mspass/seismic/CoreTimeSeries.h
+++ b/cxx/include/mspass/seismic/CoreTimeSeries.h
@@ -144,8 +144,8 @@ This method is complicated by the need to sync the changed value with
   the ends of the rhs are inside the range of the lhs.
 
   \param d is other signal to add to this.
-  \exception MsPASSError can be thrown if lhs and rhs do not have matching
-  time standards.
+  \exception MsPASSError if lhs and rhs have incompatible time standards,
+  sample intervals, or start-time grids.
   **/
   CoreTimeSeries &operator+=(const CoreTimeSeries &d);
   /*! Addition operator.
@@ -172,8 +172,8 @@ This method is complicated by the need to sync the changed value with
   the ends of the rhs are inside the range of the lhs.
 
   \param d is other signal to subract from this.
-  \exception MsPASSError can be thrown if lhs and rhs do not have matching
-  time standards.
+  \exception MsPASSError if lhs and rhs have incompatible time standards,
+  sample intervals, or start-time grids.
   **/
   CoreTimeSeries &operator-=(const CoreTimeSeries &d);
   /*! Subtraction operator.
diff --git a/cxx/src/lib/seismic/CoreSeismogram.cc b/cxx/src/lib/seismic/CoreSeismogram.cc
index e4cedee251..93a0b4253e 100644
--- a/cxx/src/lib/seismic/CoreSeismogram.cc
+++ b/cxx/src/lib/seismic/CoreSeismogram.cc
@@ -3,6 +3,7 @@
 #include "mspass/seismic/keywords.h"
 #include "mspass/utility/MsPASSError.h"
 #include "mspass/utility/SphericalCoordinate.h"
+#include "waveform_arithmetic.h"
 #include <boost/any.hpp>
 #include <float.h>
 #include <math.h>
@@ -913,39 +914,17 @@ CoreSeismogram &CoreSeismogram::operator*=(const double scale) {
   return (*this);
 }
 CoreSeismogram &CoreSeismogram::operator+=(const CoreSeismogram &d) {
-  int i, iend, jend;
-  size_t j, i0, j0;
   // Silently do nothing if d or lhs is marked dead
   if (d.dead() || (this->dead()))
     return (*this);
-  // Silently do nothing if d does not overlap with data to contain sum
-  if ((d.endtime() < mt0) || (d.mt0 > (this->endtime())))
-    return (*this);
-  if (d.tref != (this->tref))
-    throw MsPASSError("CoreSeismogram += operator cannot handle data with "
-                      "inconsistent time base\n",
-                      ErrorSeverity::Invalid);
-  /* this defines the range of left and right hand sides to be summed */
-  i = d.sample_number(this->mt0);
-  if (i < 0) {
-    j0 = this->sample_number(d.t0());
-    i0 = 0;
-  } else {
-    j0 = 0;
-    i0 = i;
-  }
-  iend = d.sample_number(this->endtime());
-  jend = this->sample_number(d.endtime());
-  if (iend >= (d.npts())) {
-    iend = d.npts() - 1;
-  }
-  if (jend >= this->npts()) {
-    jend = this->npts() - 1;
-  }
-  for (i = i0, j = j0; i <= iend && j <= jend; ++i, ++j) {
-    this->u(0, j) += d.u(0, i);
-    this->u(1, j) += d.u(1, i);
-    this->u(2, j) += d.u(2, i);
+  const auto overlap =
+      detail::arithmetic_overlap(*this, d, "CoreSeismogram::operator+=");
+  for (std::size_t i = 0; i < overlap.count; ++i) {
+    const std::size_t lhs_index = overlap.lhs_begin + i;
+    const std::size_t rhs_index = overlap.rhs_begin + i;
+    this->u(0, lhs_index) += d.u(0, rhs_index);
+    this->u(1, lhs_index) += d.u(1, rhs_index);
+    this->u(2, lhs_index) += d.u(2, rhs_index);
   }
   return (*this);
 }
@@ -953,42 +932,18 @@ CoreSeismogram &CoreSeismogram::operator+=(const CoreSeismogram &d) {
 except the += in the last loop becomes -=.  Any changes in operator+=
 must have exactly the same change here (other than a message with a
 tag to the function)*/
-CoreSeismogram &CoreSeismogram::operator-=(const CoreSeismogram &data) {
-  int i, iend, jend;
-  size_t j, i0, j0;
-  // Sun's compiler complains about const objects without this.
-  CoreSeismogram &d = const_cast<CoreSeismogram &>(data);
+CoreSeismogram &CoreSeismogram::operator-=(const CoreSeismogram &d) {
   // Silently do nothing if d is marked dead
-  if (!d.mlive)
+  if (d.dead())
     return (*this);
-  // Silently do nothing if d does not overlap with data to contain sum
-  if ((d.endtime() < mt0) || (d.mt0 > (this->endtime())))
-    return (*this);
-  if (d.tref != (this->tref))
-    throw MsPASSError("CoreSeismogram += operator cannot handle data with "
-                      "inconsistent time base\n",
-                      ErrorSeverity::Invalid);
-  /* this defines the range of left and right hand sides to be summed */
-  i = d.sample_number(this->mt0);
-  if (i < 0) {
-    j0 = this->sample_number(d.t0());
-    i0 = 0;
-  } else {
-    j0 = 0;
-    i0 = i;
-  }
-  iend = d.sample_number(this->endtime());
-  jend = this->sample_number(d.endtime());
-  if (iend >= (d.npts())) {
-    iend = d.npts() - 1;
-  }
-  if (jend >= this->npts()) {
-    jend = this->npts() - 1;
-  }
-  for (i = i0, j = j0; i <= iend && j <= jend; ++i, ++j) {
-    this->u(0, j) -= d.u(0, i);
-    this->u(1, j) -= d.u(1, i);
-    this->u(2, j) -= d.u(2, i);
+  const auto overlap =
+      detail::arithmetic_overlap(*this, d, "CoreSeismogram::operator-=");
+  for (std::size_t i = 0; i < overlap.count; ++i) {
+    const std::size_t lhs_index = overlap.lhs_begin + i;
+    const std::size_t rhs_index = overlap.rhs_begin + i;
+    this->u(0, lhs_index) -= d.u(0, rhs_index);
+    this->u(1, lhs_index) -= d.u(1, rhs_index);
+    this->u(2, lhs_index) -= d.u(2, rhs_index);
   }
   return (*this);
 }
diff --git a/cxx/src/lib/seismic/CoreTimeSeries.cc b/cxx/src/lib/seismic/CoreTimeSeries.cc
index b4924e4fa6..d557d6c8dd 100644
--- a/cxx/src/lib/seismic/CoreTimeSeries.cc
+++ b/cxx/src/lib/seismic/CoreTimeSeries.cc
@@ -3,6 +3,7 @@
 #include "mspass/seismic/keywords.h"
 #include "mspass/utility/Metadata.h"
 #include "mspass/utility/MsPASSError.h"
+#include "waveform_arithmetic.h"
 #include <vector>
 namespace mspass::seismic {
 using namespace std;
@@ -65,42 +66,13 @@ CoreTimeSeries &CoreTimeSeries::operator=(const CoreTimeSeries &tsi) {
 /*  Sum operator for CoreTimeSeries object */
 
 CoreTimeSeries &CoreTimeSeries::operator+=(const CoreTimeSeries &data) {
-  int i, iend, jend;
-  size_t i0;
-  size_t j, j0;
-  // Sun's compiler complains about const objects without this.
-  CoreTimeSeries &d = const_cast<CoreTimeSeries &>(data);
-  // Silently do nothing if d is marked dead
-  if (!d.mlive)
+  // Silently do nothing if data is marked dead
+  if (data.dead())
     return (*this);
-  // Silently do nothing if d does not overlap with data to contain sum
-  if ((d.endtime() < mt0) || (d.mt0 > (this->endtime())))
-    return (*this);
-  if (d.tref != (this->tref))
-    throw MsPASSError("CoreTimeSeries += operator cannot handle data with "
-                      "inconsistent time base\n",
-                      ErrorSeverity::Invalid);
-  /* this defines the range of left and right hand sides to be summed */
-  i = d.sample_number(this->mt0);
-  if (i < 0) {
-    j0 = this->sample_number(d.t0());
-    i0 = 0;
-  } else {
-    j0 = 0;
-    i0 = i;
-  }
-  iend = d.sample_number(this->endtime());
-  jend = this->sample_number(d.endtime());
-  if (iend >= (d.npts())) {
-    iend = d.npts() - 1;
-  }
-  if (jend >= this->npts()) {
-    jend = this->npts() - 1;
-  }
-  // cout << "i0="<<i0<<" j0="<<j0<<" iend="<<iend<<" jend="<<jend<<endl;
-  /*  Now do the actual sum using the computed ranges */
-  for (i = i0, j = j0; i <= iend && j <= jend; ++i, ++j)
-    this->s[j] += d.s[i];
+  const auto overlap =
+      detail::arithmetic_overlap(*this, data, "CoreTimeSeries::operator+=");
+  for (std::size_t i = 0; i < overlap.count; ++i)
+    this->s[overlap.lhs_begin + i] += data.s[overlap.rhs_begin + i];
   return (*this);
 }
 /* IMPORTANT:  this code is absolutely identical to that for operator+=
@@ -108,42 +80,13 @@ except the += in the last loop becomes -=.  Any changes in operator+=
 must have exactly the same change here (other than a message with a
 tag to the function)*/
 CoreTimeSeries &CoreTimeSeries::operator-=(const CoreTimeSeries &data) {
-  int i, iend, jend;
-  size_t i0;
-  size_t j, j0;
-  // Sun's compiler complains about const objects without this.
-  CoreTimeSeries &d = const_cast<CoreTimeSeries &>(data);
-  // Silently do nothing if d is marked dead
-  if (!d.mlive)
+  // Silently do nothing if data is marked dead
+  if (data.dead())
     return (*this);
-  // Silently do nothing if d does not overlap with data to contain sum
-  if ((d.endtime() < mt0) || (d.mt0 > (this->endtime())))
-    return (*this);
-  if (d.tref != (this->tref))
-    throw MsPASSError("CoreTimeSeries += operator cannot handle data with "
-                      "inconsistent time base\n",
-                      ErrorSeverity::Invalid);
-  /* this defines the range of left and right hand sides to be summed */
-  i = d.sample_number(this->mt0);
-  if (i < 0) {
-    j0 = this->sample_number(d.t0());
-    i0 = 0;
-  } else {
-    j0 = 0;
-    i0 = i;
-  }
-  iend = d.sample_number(this->endtime());
-  jend = this->sample_number(d.endtime());
-  if (iend >= (d.npts())) {
-    iend = d.npts() - 1;
-  }
-  if (jend >= this->npts()) {
-    jend = this->npts() - 1;
-  }
-  // cout << "i0="<<i0<<" j0="<<j0<<" iend="<<iend<<" jend="<<jend<<endl;
-  /*  Now do the actual sum using the computed ranges */
-  for (i = i0, j = j0; i <= iend && j <= jend; ++i, ++j)
-    this->s[j] -= d.s[i];
+  const auto overlap =
+      detail::arithmetic_overlap(*this, data, "CoreTimeSeries::operator-=");
+  for (std::size_t i = 0; i < overlap.count; ++i)
+    this->s[overlap.lhs_begin + i] -= data.s[overlap.rhs_begin + i];
   return (*this);
 }
 const CoreTimeSeries
diff --git a/cxx/src/lib/seismic/waveform_arithmetic.h b/cxx/src/lib/seismic/waveform_arithmetic.h
new file mode 100644
index 0000000000..d8e9ed788b
--- /dev/null
+++ b/cxx/src/lib/seismic/waveform_arithmetic.h
@@ -0,0 +1,67 @@
+#ifndef MSPASS_SEISMIC_WAVEFORM_ARITHMETIC_H
+#define MSPASS_SEISMIC_WAVEFORM_ARITHMETIC_H
+
+#include "mspass/seismic/BasicTimeSeries.h"
+#include "mspass/utility/MsPASSError.h"
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <sstream>
+#include <string>
+
+namespace mspass::seismic::detail {
+struct ArithmeticOverlap {
+  std::size_t lhs_begin;
+  std::size_t rhs_begin;
+  std::size_t count;
+};
+
+inline ArithmeticOverlap arithmetic_overlap(const BasicTimeSeries &lhs,
+                                            const BasicTimeSeries &rhs,
+                                            const char *caller) {
+  using mspass::utility::ErrorSeverity;
+  using mspass::utility::MsPASSError;
+
+  if (lhs.timetype() != rhs.timetype())
+    throw MsPASSError(std::string(caller) +
+                          ": operands use inconsistent time references",
+                      ErrorSeverity::Invalid);
+
+  const double lhs_dt = lhs.dt();
+  const double rhs_dt = rhs.dt();
+  const double dt_tolerance =
+      1.0e-6 * std::max(std::abs(lhs_dt), std::abs(rhs_dt));
+  if (!std::isfinite(lhs_dt) || !std::isfinite(rhs_dt) ||
+      !(std::abs(lhs_dt - rhs_dt) <= dt_tolerance)) {
+    std::ostringstream message;
+    message << caller << ": sample intervals do not match: lhs dt=" << lhs_dt
+            << ", rhs dt=" << rhs_dt;
+    throw MsPASSError(message.str(), ErrorSeverity::Invalid);
+  }
+
+  const double offset_samples = (rhs.t0() - lhs.t0()) / lhs_dt;
+  const double rounded_offset = std::round(offset_samples);
+  if (!(std::abs(offset_samples - rounded_offset) <= 1.0e-6)) {
+    std::ostringstream message;
+    message << caller << ": start times are not aligned to the sample grid: "
+            << "offset=" << offset_samples << " samples";
+    throw MsPASSError(message.str(), ErrorSeverity::Invalid);
+  }
+
+  if (rounded_offset >= static_cast<double>(lhs.npts()) ||
+      rounded_offset <= -static_cast<double>(rhs.npts()))
+    return ArithmeticOverlap{0, 0, 0};
+
+  if (rounded_offset >= 0.0) {
+    const std::size_t lhs_begin = static_cast<std::size_t>(rounded_offset);
+    return ArithmeticOverlap{lhs_begin, 0,
+                             std::min(lhs.npts() - lhs_begin, rhs.npts())};
+  }
+
+  const std::size_t rhs_begin = static_cast<std::size_t>(-rounded_offset);
+  return ArithmeticOverlap{0, rhs_begin,
+                           std::min(lhs.npts(), rhs.npts() - rhs_begin)};
+}
+} // namespace mspass::seismic::detail
+
+#endif
diff --git a/cxx/test/CMakeLists.txt b/cxx/test/CMakeLists.txt
index 4a3b795346..e972974af7 100644
--- a/cxx/test/CMakeLists.txt
+++ b/cxx/test/CMakeLists.txt
@@ -1,6 +1,7 @@
 include(CTest)
 if(BUILD_TESTING)
   add_subdirectory(amap)
+  add_subdirectory(arithmetic)
   add_subdirectory(decon_serialization)
   add_subdirectory(data_directory)
   add_subdirectory(dmatrix)
@@ -19,6 +20,8 @@ if(BUILD_TESTING)
 
 # note this test needs to access a data file in the run directory. 
 # we set that with WORKING_DIRECTORY and the special symbol after the keyword
+  add_test(NAME test_waveform_arithmetic
+    COMMAND ${PROJECT_BINARY_DIR}/test/arithmetic/test_waveform_arithmetic)
   add_test(NAME test_decon_serialization COMMAND ${PROJECT_BINARY_DIR}/test/decon_serialization/test_decon_serialization WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/decon_serialization)
   add_test(NAME test_data_directory
     COMMAND ${PROJECT_BINARY_DIR}/test/data_directory/test_data_directory)
diff --git a/cxx/test/arithmetic/CMakeLists.txt b/cxx/test/arithmetic/CMakeLists.txt
new file mode 100644
index 0000000000..1cd4b5659a
--- /dev/null
+++ b/cxx/test/arithmetic/CMakeLists.txt
@@ -0,0 +1,8 @@
+add_executable(test_waveform_arithmetic test_waveform_arithmetic.cc)
+include_directories(
+  ${Boost_INCLUDE_DIRS}
+  ${pybind11_INCLUDE_DIR}
+  ${Python_INCLUDE_DIRS}
+  ${PROJECT_SOURCE_DIR}/include/)
+
+target_link_libraries(test_waveform_arithmetic PRIVATE mspass ${Boost_LIBRARIES} ${Python_LIBRARIES})
diff --git a/cxx/test/arithmetic/test_waveform_arithmetic.cc b/cxx/test/arithmetic/test_waveform_arithmetic.cc
new file mode 100644
index 0000000000..3c565dcb5c
--- /dev/null
+++ b/cxx/test/arithmetic/test_waveform_arithmetic.cc
@@ -0,0 +1,235 @@
+#include "mspass/seismic/Seismogram.h"
+#include "mspass/seismic/TimeSeries.h"
+#include "mspass/utility/MsPASSError.h"
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <iostream>
+#include <limits>
+#include <string>
+#include <vector>
+
+using mspass::seismic::Seismogram;
+using mspass::seismic::TimeReferenceType;
+using mspass::seismic::TimeSeries;
+using mspass::utility::ErrorSeverity;
+using mspass::utility::MsPASSError;
+
+namespace {
+int failure_count{0};
+
+void check(const bool condition, const std::string &message) {
+  if (!condition) {
+    std::cerr << message << '\n';
+    ++failure_count;
+  }
+}
+
+template <class Waveform> struct Samples;
+
+template <> struct Samples<TimeSeries> {
+  static constexpr std::size_t component_count{1};
+  static double get(const TimeSeries &data, const std::size_t,
+                    const std::size_t sample) {
+    return data.s[sample];
+  }
+  static void set(TimeSeries &data, const std::size_t, const std::size_t sample,
+                  const double value) {
+    data.s[sample] = value;
+  }
+};
+
+template <> struct Samples<Seismogram> {
+  static constexpr std::size_t component_count{3};
+  static double get(const Seismogram &data, const std::size_t component,
+                    const std::size_t sample) {
+    return data.u(component, sample);
+  }
+  static void set(Seismogram &data, const std::size_t component,
+                  const std::size_t sample, const double value) {
+    data.u(component, sample) = value;
+  }
+};
+
+template <class Waveform>
+Waveform make_waveform(const std::size_t sample_count, const double t0,
+                       const double dt, const double base) {
+  Waveform data(sample_count);
+  data.set_t0(t0);
+  data.set_dt(dt);
+  data.set_tref(TimeReferenceType::Relative);
+  data.set_live();
+  data.put("sentinel", std::string("unchanged"));
+  for (std::size_t component = 0;
+       component < Samples<Waveform>::component_count; ++component)
+    for (std::size_t sample = 0; sample < sample_count; ++sample)
+      Samples<Waveform>::set(data, component, sample,
+                             base + 10.0 * static_cast<double>(component) +
+                                 static_cast<double>(sample));
+  return data;
+}
+
+template <class Waveform> struct Snapshot {
+  std::size_t sample_count;
+  double t0;
+  double dt;
+  TimeReferenceType time_reference;
+  bool live;
+  int error_count;
+  std::string sentinel;
+  std::vector<double> samples;
+};
+
+template <class Waveform> Snapshot<Waveform> snapshot(const Waveform &data) {
+  Snapshot<Waveform> result{data.npts(),
+                            data.t0(),
+                            data.dt(),
+                            data.timetype(),
+                            data.live(),
+                            data.elog.size(),
+                            data.get_string("sentinel"),
+                            {}};
+  result.samples.reserve(Samples<Waveform>::component_count * data.npts());
+  for (std::size_t component = 0;
+       component < Samples<Waveform>::component_count; ++component)
+    for (std::size_t sample = 0; sample < data.npts(); ++sample)
+      result.samples.push_back(Samples<Waveform>::get(data, component, sample));
+  return result;
+}
+
+template <class Waveform>
+void check_state(const Waveform &actual, const Snapshot<Waveform> &expected,
+                 const std::string &context) {
+  check(actual.npts() == expected.sample_count, context + ": npts changed");
+  check(actual.t0() == expected.t0, context + ": t0 changed");
+  check(actual.dt() == expected.dt, context + ": dt changed");
+  check(actual.timetype() == expected.time_reference,
+        context + ": time reference changed");
+  check(actual.live() == expected.live, context + ": live state changed");
+  check(actual.elog.size() == expected.error_count,
+        context + ": error log changed");
+  check(actual.get_string("sentinel") == expected.sentinel,
+        context + ": metadata changed");
+  std::size_t index{0};
+  for (std::size_t component = 0;
+       component < Samples<Waveform>::component_count; ++component)
+    for (std::size_t sample = 0; sample < actual.npts(); ++sample, ++index)
+      check(Samples<Waveform>::get(actual, component, sample) ==
+                expected.samples[index],
+            context + ": sample changed at component " +
+                std::to_string(component) + ", index " +
+                std::to_string(sample));
+}
+
+template <class Waveform>
+void combine(Waveform &lhs, const Waveform &rhs, const bool add) {
+  if (add)
+    lhs += rhs;
+  else
+    lhs -= rhs;
+}
+
+template <class Waveform>
+void verify_valid(const bool add, const double rhs_t0, const double rhs_dt,
+                  const int expected_offset, const std::string &context,
+                  const double lhs_dt = 1.0) {
+  Waveform lhs = make_waveform<Waveform>(5, 0.0, lhs_dt, 100.0);
+  const Waveform rhs = make_waveform<Waveform>(4, rhs_t0, rhs_dt, 10.0);
+  Waveform expected(lhs);
+  for (std::size_t rhs_index = 0; rhs_index < rhs.npts(); ++rhs_index) {
+    const int lhs_index = expected_offset + static_cast<int>(rhs_index);
+    if (lhs_index < 0 || lhs_index >= static_cast<int>(lhs.npts()))
+      continue;
+    for (std::size_t component = 0;
+         component < Samples<Waveform>::component_count; ++component) {
+      const double value = Samples<Waveform>::get(lhs, component, lhs_index) +
+                           (add ? 1.0 : -1.0) * Samples<Waveform>::get(
+                                                    rhs, component, rhs_index);
+      Samples<Waveform>::set(expected, component, lhs_index, value);
+    }
+  }
+
+  combine(lhs, rhs, add);
+  check_state(lhs, snapshot(expected), context);
+}
+
+template <class Waveform>
+void verify_rejected(const bool add, const double rhs_t0, const double rhs_dt,
+                     const std::string &context, const double lhs_dt = 1.0) {
+  Waveform lhs = make_waveform<Waveform>(5, 0.0, lhs_dt, 100.0);
+  Waveform rhs = make_waveform<Waveform>(4, rhs_t0, rhs_dt, 10.0);
+  rhs.elog.log_error("rhs", "must not merge on rejection",
+                     ErrorSeverity::Complaint);
+  const auto before = snapshot(lhs);
+  bool threw_invalid{false};
+  try {
+    combine(lhs, rhs, add);
+  } catch (const MsPASSError &error) {
+    threw_invalid = error.severity() == ErrorSeverity::Invalid;
+  }
+  check(threw_invalid, context + ": did not throw MsPASSError Invalid");
+  check_state(lhs, before, context + " rejection atomicity");
+}
+
+template <class Waveform>
+void run_suite(const bool add, const std::string &type_name) {
+  const std::string operation = add ? "+=" : "-=";
+  const std::string prefix = type_name + " " + operation;
+
+  verify_valid<Waveform>(add, 0.0, 1.0, 0, prefix + " equal grid");
+  verify_valid<Waveform>(add, 2.0, 1.0, 2, prefix + " positive offset");
+  verify_valid<Waveform>(add, -2.0, 1.0, -2, prefix + " negative offset");
+  verify_valid<Waveform>(add, 5.0, 1.0, 5, prefix + " positive no overlap");
+  verify_valid<Waveform>(add, -4.0, 1.0, -4, prefix + " negative no overlap");
+
+  constexpr double tolerance = 1.0e-6;
+  const double positive_offset_at_tolerance = 1.0 + tolerance;
+  const double positive_offset_beyond = std::nextafter(
+      positive_offset_at_tolerance, std::numeric_limits<double>::infinity());
+  const double negative_offset_at_tolerance = -1.0 - tolerance;
+  const double negative_offset_beyond = std::nextafter(
+      negative_offset_at_tolerance, -std::numeric_limits<double>::infinity());
+  check(std::abs(positive_offset_at_tolerance -
+                 std::round(positive_offset_at_tolerance)) <= tolerance,
+        "positive offset boundary fixture is invalid");
+  check(std::abs(positive_offset_beyond - std::round(positive_offset_beyond)) >
+            tolerance,
+        "positive offset rejection fixture is invalid");
+  check(std::abs(negative_offset_at_tolerance -
+                 std::round(negative_offset_at_tolerance)) <= tolerance,
+        "negative offset boundary fixture is invalid");
+  check(std::abs(negative_offset_beyond - std::round(negative_offset_beyond)) >
+            tolerance,
+        "negative offset rejection fixture is invalid");
+  verify_valid<Waveform>(add, positive_offset_at_tolerance, 1.0, 1,
+                         prefix + " positive offset at tolerance");
+  verify_valid<Waveform>(add, negative_offset_at_tolerance, 1.0, -1,
+                         prefix + " negative offset at tolerance");
+  verify_rejected<Waveform>(add, positive_offset_beyond, 1.0,
+                            prefix + " positive offset beyond tolerance");
+  verify_rejected<Waveform>(add, negative_offset_beyond, 1.0,
+                            prefix + " negative offset beyond tolerance");
+
+  constexpr double lhs_dt = 1.0e6;
+  constexpr double dt_at_tolerance = lhs_dt - 1.0;
+  const double dt_beyond = std::nextafter(dt_at_tolerance, 0.0);
+  check(std::abs(lhs_dt - dt_at_tolerance) ==
+            tolerance * std::max(std::abs(lhs_dt), std::abs(dt_at_tolerance)),
+        "dt boundary fixture is invalid");
+  check(std::abs(lhs_dt - dt_beyond) >
+            tolerance * std::max(std::abs(lhs_dt), std::abs(dt_beyond)),
+        "dt rejection fixture is invalid");
+  verify_valid<Waveform>(add, 0.0, dt_at_tolerance, 0,
+                         prefix + " dt at relative tolerance", lhs_dt);
+  verify_rejected<Waveform>(add, 0.0, dt_beyond,
+                            prefix + " dt beyond relative tolerance", lhs_dt);
+}
+} // namespace
+
+int main() {
+  run_suite<TimeSeries>(true, "TimeSeries");
+  run_suite<TimeSeries>(false, "TimeSeries");
+  run_suite<Seismogram>(true, "Seismogram");
+  run_suite<Seismogram>(false, "Seismogram");
+  return failure_count == 0 ? 0 : 1;
+}
diff --git a/python/mspasspy/algorithms/MCXcorStacking.py b/python/mspasspy/algorithms/MCXcorStacking.py
index d526215b02..11c0e9d7df 100644
--- a/python/mspasspy/algorithms/MCXcorStacking.py
+++ b/python/mspasspy/algorithms/MCXcorStacking.py
@@ -1084,8 +1084,10 @@ def beam_align(ensemble, beam, window=None, time_shift_limit=10.0):
         because the value can to sent to a C++ method that it type
         sensitive)
     :return: the input ensemble with its members time shifted in place to align
-        with the time base of beam.  If a window is defined it is used only for
-        correlation and is not applied to the ensemble members.
+        with the time base of beam.  Each shifted start time is placed on the
+        nearest sample of the beam grid without changing any sample values.
+        If a window is defined it is used only for correlation and is not
+        applied to the ensemble members.
 
     """
     # this may not be necessary for internal use but if used
@@ -1126,6 +1128,7 @@ def beam_align(ensemble, beam, window=None, time_shift_limit=10.0):
             # not shifted from UTC times
             # also note a +lag requires a - shift
             ensemble.member[i].shift(timelag)
+            _snap_start_time_to_grid(ensemble.member[i], beam)
     return ensemble
 
 
@@ -1645,6 +1648,7 @@ def align_and_stack(
                 # this method alters the t0 values of the ensemble members
                 # when used plots will tend to be aligned with 0 relative time
                 ensemble.member[i].shift(tshift_mean)
+                _snap_start_time_to_grid(ensemble.member[i], beam)
     for i in range(len(ensemble.member)):
         if ensemble.member[i].live:
             # in this context it0_key should always be defined
@@ -2452,6 +2456,17 @@ def _xcor_shift(ts, beam):
     return lagtime
 
 
+def _snap_start_time_to_grid(datum, reference):
+    """Move ``datum.t0`` to the nearest sample on ``reference``'s grid.
+
+    This is an explicit MCXcor alignment operation:  it changes only the time
+    assigned to sample zero and never changes the sample values, ``dt``, or
+    ``npts``.  General waveform arithmetic remains strict and rejects operands
+    that have not first been placed on a common grid.
+    """
+    datum.t0 = reference.time(reference.sample_number(datum.t0))
+
+
 def _update_xcor_beam(xcorens, beam0, robust_stack_method, wts) -> TimeSeries:
     """
     Internal method used to update the beam signal used for cross correlation
@@ -2503,7 +2518,7 @@ def _update_xcor_beam(xcorens, beam0, robust_stack_method, wts) -> TimeSeries:
                     xcorens.member[i], stime, etime, short_segment_handling="pad"
                 )
                 d *= wt
-                # TimeSeries::operator+= handles time correctly so indexing is not needed here
+                _snap_start_time_to_grid(d, beam)
                 beam += d
                 sumwt += wt
         if sumwt > 0.0:
diff --git a/python/tests/algorithms/test_MCXcorStacking.py b/python/tests/algorithms/test_MCXcorStacking.py
index b8e412d230..6864271a9b 100755
--- a/python/tests/algorithms/test_MCXcorStacking.py
+++ b/python/tests/algorithms/test_MCXcorStacking.py
@@ -40,6 +40,7 @@
     _set_phases,
     regularize_sampling,
     ensemble_time_range,
+    _snap_start_time_to_grid,
 )
 
 
@@ -391,6 +392,27 @@ def test_align_and_stack():
     validate_lags(eo, lag_in_samples)
 
 
+def test_mcxcor_nearest_grid_alignment_preserves_samples():
+    reference = TimeSeries(8)
+    reference.dt = 0.05
+    reference.t0 = -1.0
+    reference.set_live()
+
+    datum = TimeSeries(6)
+    datum.dt = reference.dt
+    datum.t0 = reference.t0 + 2.49 * reference.dt
+    datum.set_live()
+    datum.data = DoubleVector(np.arange(datum.npts, dtype=float))
+    samples_before = np.array(datum.data, copy=True)
+
+    _snap_start_time_to_grid(datum, reference)
+
+    assert datum.t0 == reference.time(2)
+    assert datum.dt == reference.dt
+    assert datum.npts == len(samples_before)
+    np.testing.assert_array_equal(datum.data, samples_before)
+
+
 def test_align_and_stack_error_handlers():
     """
     Like above but only tests error handlers.  Uses the
diff --git a/python/tests/global_history/test_manager_no_spark_dask.py b/python/tests/global_history/test_manager_no_spark_dask.py
index 40033f959e..d931377520 100644
--- a/python/tests/global_history/test_manager_no_spark_dask.py
+++ b/python/tests/global_history/test_manager_no_spark_dask.py
@@ -2,7 +2,7 @@
 import sys
 
 sys.path.append("python/tests")
-from helper import get_live_timeseries
+from helper import get_live_timeseries, get_live_timeseries_list
 
 with mock.patch.dict(sys.modules, {"pyspark": None, "dask": None}):
     from mspasspy.reduce import stack
@@ -39,7 +39,7 @@ def test_normal_map():
             assert test_map_res[i].data == t[i].data + t[i].data
 
     def test_normal_reduce():
-        t = [get_live_timeseries() for i in range(5)]
+        t = get_live_timeseries_list(5)
         s = t[0]
         for i in range(1, 5):
             s += t[i]
diff --git a/python/tests/global_history/test_manager_spark_dask.py b/python/tests/global_history/test_manager_spark_dask.py
index 195d56654b..8a501f470f 100644
--- a/python/tests/global_history/test_manager_spark_dask.py
+++ b/python/tests/global_history/test_manager_spark_dask.py
@@ -29,6 +29,7 @@
 from helper import (
     get_live_seismogram,
     get_live_timeseries,
+    get_live_timeseries_list,
     get_live_timeseries_ensemble,
     get_live_seismogram_ensemble,
 )
@@ -161,7 +162,7 @@ def test_normal_map(self):
             assert test_map_res[i].data == t[i].data + t[i].data
 
     def test_normal_reduce(self):
-        t = [get_live_timeseries() for i in range(5)]
+        t = get_live_timeseries_list(5)
         s = t[0]
         for i in range(1, 5):
             s += t[i]
@@ -501,7 +502,7 @@ def test_mspass_reduce(self, spark_context):
         manager_db = Database(self.client, "test_manager")
         manager_db["history_global"].delete_many({})
 
-        l = [get_live_timeseries() for i in range(5)]
+        l = get_live_timeseries_list(5)
         # test mspass_reduce for spark
         spark_res = spark_reduce(l, self.manager, spark_context)
         assert (
diff --git a/python/tests/helper.py b/python/tests/helper.py
index 3ac65afc29..cf7eab8b9e 100644
--- a/python/tests/helper.py
+++ b/python/tests/helper.py
@@ -37,7 +37,7 @@ def get_live_seismogram(ts_size=255, sampling_rate=20.0):
     return seis
 
 
-def get_live_timeseries(ts_size=255, sampling_rate=20.0):
+def get_live_timeseries(ts_size=255, sampling_rate=20.0, start_time=None):
     ts = TimeSeries()
     ts.set_live()
     ts.dt = 1 / sampling_rate
@@ -46,7 +46,9 @@ def get_live_timeseries(ts_size=255, sampling_rate=20.0):
     ts.put("npts", ts_size)
     ts.put("sampling_rate", sampling_rate)
     ts.tref = TimeReferenceType.UTC
-    ts.t0 = datetime.datetime.now().timestamp()
+    if start_time is None:
+        start_time = datetime.datetime.now().timestamp()
+    ts.t0 = start_time
     ts["delta"] = 0.1
     ts["calib"] = 0.1
     ts["site_id"] = bson.objectid.ObjectId()
@@ -90,19 +92,23 @@ def get_live_seismogram_ensemble(n):
     return seis_e
 
 
-def get_live_timeseries_ensemble(n):
+def get_live_timeseries_ensemble(n, start_time=None):
     tse = TimeSeriesEnsemble()
+    if start_time is None:
+        start_time = datetime.datetime.now().timestamp()
     for i in range(n):
-        ts = get_live_timeseries()
+        ts = get_live_timeseries(start_time=start_time)
         tse.member.append(ts)
     tse.set_live()
     return tse
 
 
-def get_live_timeseries_list(n, ts_size=255):
+def get_live_timeseries_list(n, ts_size=255, start_time=None):
     ts_e = []
+    if start_time is None:
+        start_time = datetime.datetime.now().timestamp()
     for i in range(n):
-        ts = get_live_timeseries(ts_size)
+        ts = get_live_timeseries(ts_size, start_time=start_time)
         ts_e.append(ts)
     return ts_e
 
diff --git a/python/tests/test_reduce.py b/python/tests/test_reduce.py
index c40a7f4c37..6f7f47eb3f 100644
--- a/python/tests/test_reduce.py
+++ b/python/tests/test_reduce.py
@@ -25,6 +25,7 @@
 from helper import (
     get_live_seismogram,
     get_live_timeseries,
+    get_live_timeseries_list,
     get_live_timeseries_ensemble,
     get_live_seismogram_ensemble,
     get_live_seismogram_list,
@@ -84,14 +85,15 @@ def test_reduce_stack():
     for i in range(3):
         assert np.isclose(seis1.data[i], res[i]).all()  # fixme
 
-    ts1 = get_live_timeseries()
-    ts2 = get_live_timeseries()
+    start_time = 1234567890.0
+    ts1 = get_live_timeseries(start_time=start_time)
+    ts2 = get_live_timeseries(start_time=start_time)
     ts1_cp = np.array(ts1.data)
     stack(ts1, ts2)
     assert np.isclose(ts1.data, (np.array(ts1_cp) + np.array(ts2.data))).all()
 
-    tse1 = get_live_timeseries_ensemble(2)
-    tse2 = get_live_timeseries_ensemble(2)
+    tse1 = get_live_timeseries_ensemble(2, start_time=start_time)
+    tse2 = get_live_timeseries_ensemble(2, start_time=start_time)
     tse1_cp = TimeSeriesEnsemble(tse1)
     stack(tse1, tse2)
     for i in range(2):
@@ -112,6 +114,23 @@ def test_reduce_stack():
             assert np.isclose(seis_e1.member[i].data[j], res[j]).all()  # fixme
 
 
+def test_timeseries_helpers_share_explicit_start_time():
+    start_time = 0.0
+    timeseries_list = get_live_timeseries_list(3, start_time=start_time)
+    timeseries_ensemble = get_live_timeseries_ensemble(3, start_time=start_time)
+
+    assert [datum.t0 for datum in timeseries_list] == [start_time] * 3
+    assert [datum.t0 for datum in timeseries_ensemble.member] == [start_time] * 3
+    list_other_values = [datum.data[0] for datum in timeseries_list[1:]]
+    ensemble_other_values = [datum.data[0] for datum in timeseries_ensemble.member[1:]]
+    timeseries_list[0].data[0] = 1.0
+    timeseries_ensemble.member[0].data[0] = 1.0
+    assert [datum.data[0] for datum in timeseries_list[1:]] == list_other_values
+    assert [
+        datum.data[0] for datum in timeseries_ensemble.member[1:]
+    ] == ensemble_other_values
+
+
 def test_reduce_stack_exception():
     tse1 = get_live_timeseries_ensemble(2)
     tse2 = get_live_timeseries_ensemble(3)
@@ -144,7 +163,7 @@ def spark_reduce(input, sc):
 
 
 def test_reduce_dask_spark(spark_context):
-    l = [get_live_timeseries() for i in range(5)]
+    l = get_live_timeseries_list(5)
     res = np.zeros(255)
     for i in range(5):
         for j in range(255):
diff --git a/python/tests/test_waveform_arithmetic.py b/python/tests/test_waveform_arithmetic.py
new file mode 100644
index 0000000000..611014a5b9
--- /dev/null
+++ b/python/tests/test_waveform_arithmetic.py
@@ -0,0 +1,139 @@
+import math
+
+import numpy as np
+import pytest
+
+from mspasspy.ccore.seismic import Seismogram, TimeReferenceType, TimeSeries
+from mspasspy.ccore.utility import ErrorSeverity, MsPASSError
+
+
+def _make_waveform(waveform_type, sample_count, t0, dt, base):
+    waveform = waveform_type(sample_count)
+    waveform.t0 = t0
+    waveform.dt = dt
+    waveform.tref = TimeReferenceType.Relative
+    waveform.set_live()
+    waveform["sentinel"] = "unchanged"
+    if waveform_type is TimeSeries:
+        for sample in range(sample_count):
+            waveform.data[sample] = base + sample
+    else:
+        for component in range(3):
+            waveform.data[component, :] = (
+                base + 10.0 * component + np.arange(sample_count)
+            )
+    return waveform
+
+
+def _snapshot(waveform):
+    return {
+        "npts": waveform.npts,
+        "t0": waveform.t0,
+        "dt": waveform.dt,
+        "tref": waveform.tref,
+        "live": waveform.live,
+        "error_count": waveform.elog.size(),
+        "sentinel": waveform["sentinel"],
+        "data": np.array(waveform.data, copy=True),
+    }
+
+
+def _assert_state(waveform, expected):
+    assert waveform.npts == expected["npts"]
+    assert waveform.t0 == expected["t0"]
+    assert waveform.dt == expected["dt"]
+    assert waveform.tref == expected["tref"]
+    assert waveform.live == expected["live"]
+    assert waveform.elog.size() == expected["error_count"]
+    assert waveform["sentinel"] == expected["sentinel"]
+    np.testing.assert_array_equal(waveform.data, expected["data"])
+
+
+def _combine(lhs, rhs, operation):
+    if operation == "add":
+        lhs += rhs
+    else:
+        lhs -= rhs
+
+
+def _verify_valid(waveform_type, operation, rhs_t0, rhs_dt, offset, lhs_dt=1.0):
+    lhs = _make_waveform(waveform_type, 5, 0.0, lhs_dt, 100.0)
+    rhs = _make_waveform(waveform_type, 4, rhs_t0, rhs_dt, 10.0)
+    expected = _snapshot(lhs)
+    sign = 1.0 if operation == "add" else -1.0
+    for rhs_index in range(rhs.npts):
+        lhs_index = offset + rhs_index
+        if 0 <= lhs_index < lhs.npts:
+            if waveform_type is TimeSeries:
+                expected["data"][lhs_index] += sign * rhs.data[rhs_index]
+            else:
+                expected["data"][:, lhs_index] += sign * rhs.data[:, rhs_index]
+
+    _combine(lhs, rhs, operation)
+    _assert_state(lhs, expected)
+
+
+def _verify_rejected(waveform_type, operation, rhs_t0, rhs_dt, lhs_dt=1.0):
+    lhs = _make_waveform(waveform_type, 5, 0.0, lhs_dt, 100.0)
+    rhs = _make_waveform(waveform_type, 4, rhs_t0, rhs_dt, 10.0)
+    rhs.elog.log_error("rhs", "must not merge on rejection", ErrorSeverity.Complaint)
+    before = _snapshot(lhs)
+
+    with pytest.raises(MsPASSError) as exc_info:
+        _combine(lhs, rhs, operation)
+
+    assert exc_info.value.severity == ErrorSeverity.Invalid
+    _assert_state(lhs, before)
+
+
+@pytest.mark.parametrize("waveform_type", [TimeSeries, Seismogram])
+@pytest.mark.parametrize("operation", ["add", "subtract"])
+def test_waveform_arithmetic_grid_contract(waveform_type, operation):
+    valid_cases = [
+        (0.0, 1.0, 0),
+        (2.0, 1.0, 2),
+        (-2.0, 1.0, -2),
+        (5.0, 1.0, 5),
+        (-4.0, 1.0, -4),
+    ]
+
+    tolerance = 1.0e-6
+    positive_offset_at_tolerance = 1.0 + tolerance
+    negative_offset_at_tolerance = -1.0 - tolerance
+    assert (
+        abs(positive_offset_at_tolerance - round(positive_offset_at_tolerance))
+        <= tolerance
+    )
+    assert (
+        abs(negative_offset_at_tolerance - round(negative_offset_at_tolerance))
+        <= tolerance
+    )
+    valid_cases.extend(
+        [
+            (positive_offset_at_tolerance, 1.0, 1),
+            (negative_offset_at_tolerance, 1.0, -1),
+        ]
+    )
+
+    for rhs_t0, rhs_dt, offset in valid_cases:
+        _verify_valid(waveform_type, operation, rhs_t0, rhs_dt, offset)
+
+    lhs_dt = 1.0e6
+    dt_at_tolerance = lhs_dt - 1.0
+    dt_beyond = math.nextafter(dt_at_tolerance, 0.0)
+    assert abs(lhs_dt - dt_at_tolerance) == tolerance * max(
+        abs(lhs_dt), abs(dt_at_tolerance)
+    )
+    assert abs(lhs_dt - dt_beyond) > tolerance * max(abs(lhs_dt), abs(dt_beyond))
+    _verify_valid(waveform_type, operation, 0.0, dt_at_tolerance, 0, lhs_dt=lhs_dt)
+
+    positive_offset_beyond = math.nextafter(positive_offset_at_tolerance, math.inf)
+    negative_offset_beyond = math.nextafter(negative_offset_at_tolerance, -math.inf)
+    assert abs(positive_offset_beyond - round(positive_offset_beyond)) > tolerance
+    assert abs(negative_offset_beyond - round(negative_offset_beyond)) > tolerance
+    for rhs_t0, rhs_dt in [
+        (positive_offset_beyond, 1.0),
+        (negative_offset_beyond, 1.0),
+    ]:
+        _verify_rejected(waveform_type, operation, rhs_t0, rhs_dt)
+    _verify_rejected(waveform_type, operation, 0.0, dt_beyond, lhs_dt=lhs_dt)