diff --git a/HISTORY.rst b/HISTORY.rst
index b07eb83..3ec7265 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -2,12 +2,14 @@
 History
 =======
 
-1.1.X (2019-10-XX)
+1.1.1 (2019-10-29)
 -------------------
 
 * Fixed issue in `get_zero_crossings_array_indices` where it would fail if array did not contain any zeros.
 * Added calculation of equivalent number of cycles and equivalent uniform amplitude using power law relationship as intensity measures
 * Added function `get_n_cyc_array` to compute number of cycles series from a loading series
+* Added intensity measure `im.calc_unit_kinetic_energy()` to compute the cumulative change in kinetic energy according to Millen et al. (2019)
+* Added `surface.py` with calculation of surface energy and cumulative change in surface energy time series versus depth from surface
 
 
 1.1.0 (2019-10-08)
diff --git a/eqsig/__about__.py b/eqsig/__about__.py
index 36c0ef7..bc7f771 100644
--- a/eqsig/__about__.py
+++ b/eqsig/__about__.py
@@ -1,4 +1,4 @@
 __project__ = "eqsig"
 __author__ = "Maxim Millen"
-__version__ = "1.1.0"
+__version__ = "1.1.1"
 __license__ = "MIT"
diff --git a/eqsig/__init__.py b/eqsig/__init__.py
index 6b55422..10b501c 100644
--- a/eqsig/__init__.py
+++ b/eqsig/__init__.py
@@ -10,6 +10,7 @@
 from eqsig import loader
 from eqsig import __about__
 from eqsig import sdof
+from eqsig import surface
 
 __project__ = __about__.__project__
 __author__ = __about__.__author__
diff --git a/eqsig/im.py b/eqsig/im.py
index 4d9acae..6845e3b 100644
--- a/eqsig/im.py
+++ b/eqsig/im.py
@@ -524,3 +524,22 @@ def calc_cyc_amp_combined_arrays_w_power_law(values0, values1, n_cyc, b):
     csr_n15_series = np.cumsum((np.abs(csr_peaks_s0) ** (1. / b) + np.abs(csr_peaks_s1) ** (1. / b)) / 2 / n_cyc) ** b
 
     return csr_n15_series
+
+
+def calc_unit_kinetic_energy(acc_signal):
+    """
+    Calculates the cumulative absolute change in kinetic energy for a unit volume of soil with unit mass
+
+    Parameters
+    ----------
+    acc_signal: eqsig.AccSignal
+
+    Returns
+    -------
+
+    """
+    kin_energy = 0.5 * acc_signal.velocity * np.abs(acc_signal.velocity)
+    delta_energy = np.diff(kin_energy)
+    delta_energy = np.insert(delta_energy, 0, kin_energy[0])
+    cum_delta_energy = np.cumsum(abs(delta_energy))
+    return cum_delta_energy
diff --git a/eqsig/sdof.py b/eqsig/sdof.py
index 7603b79..a7d7138 100644
--- a/eqsig/sdof.py
+++ b/eqsig/sdof.py
@@ -223,6 +223,33 @@ def time_the_generation_of_response_spectra():
     s_d, s_v, s_a = pseudo_response_spectra(motion, step, periods, xi)
 
 
+def calc_resp_uke_spectrum(acc_signal, periods=None, xi=None):
+    """
+    Calculates the sdof response (kinematic + stored) energy spectrum
+
+    :param acc_signal:
+    :param periods:
+    :param xi:
+    :return:
+    """
+
+    if periods is None:
+        periods = acc_signal.response_times
+    else:
+        periods = np.array(periods)
+    if xi is None:
+        xi = 0.05
+    resp_u, resp_v, resp_a = response_series(acc_signal.values, acc_signal.dt, periods, xi)
+    mass = 1
+
+    kin_energy = 0.5 * resp_v ** 2 * mass
+    delta_energy = np.diff(kin_energy)
+    # double for strain then half since only positive increasing
+    cum_delta_energy = np.sum(abs(delta_energy), axis=1)
+
+    return cum_delta_energy
+
+
 if __name__ == '__main__':
 
     # time_response_spectra()
diff --git a/eqsig/surface.py b/eqsig/surface.py
new file mode 100644
index 0000000..a98a156
--- /dev/null
+++ b/eqsig/surface.py
@@ -0,0 +1,103 @@
+import numpy as np
+import scipy.integrate
+from eqsig import functions as fn
+
+
+def trim_to_length(values, npts, surf2depth_travel_times, dt, trim=False, start=False, s2s_travel_time=0.0):
+    """
+    Trim a 2D array back to a have the same length
+
+    Parameters
+    ----------
+    values: array_like
+        Values to be outputted
+    s2s_travel_time: float
+        travel time from input motion location to the surface
+    surf2depth_travel_times: array_like
+        travel time from surface to depths of interest
+    trim: bool
+        if true then forces array to be same length as values array
+    start: bool
+        if true then forces array to have the same start time as values array
+
+    Returns
+    -------
+
+    """
+
+    surf_to_depth_shifts = np.array(surf2depth_travel_times / dt, dtype=int)
+    start_shift = int(s2s_travel_time / dt)
+    if start:  # Trim front
+        sis = start_shift - surf_to_depth_shifts
+        if not trim:
+            extras = np.max([sis, 0]) - np.min([np.min(2 * surf_to_depth_shifts), 0])
+            npts = npts + extras  # plus the zero padding in front and back
+    else:
+        if trim:
+            sis = np.zeros_like(surf_to_depth_shifts)
+        else:  # no changes required
+            return values
+    outs = np.zeros((len(surf_to_depth_shifts), npts))
+    print('len_outs: ', outs.shape)
+    for i in range(len(surf_to_depth_shifts)):
+        if sis[i] < 0:
+            outs[i] = values[i, -sis[i]: npts - sis[i]]
+        else:
+            outs[i, sis[i]:] = values[i, : npts - sis[i]]  # zero padded
+    return outs
+
+
+def calc_surface_energy(asig, travel_times, nodal=True, up_red=1, down_red=1, stt=0.0, trim=False, start=False):
+    """
+    Calculates the energy at different travel times from a surface
+
+    Parameters
+    ----------
+    asig
+    travel_times: array_like
+        Travel times from surface to depth of interest
+    nodal: bool
+        If true then surface is nodal (minima)
+    up_red: float or array_like
+        upward wave reduction factors
+    down_red: float or array_like
+        downward wave reduction factors
+    trim: bool
+        if true then forces array to be same length as values array
+    start: bool
+        if true then forces array to have the same start time as values array
+    stt: float
+        Travel time from input motion location to the surface
+
+    Returns
+    -------
+
+    """
+    shifts = np.array(2 * travel_times / asig.dt, dtype=int)
+    up_wave = np.pad(asig.values, (0, np.max(shifts)), mode='constant', constant_values=0) * up_red  # 1d
+    down_waves = fn.put_array_in_2d_array(asig.values, shifts) * down_red
+    if nodal:
+        acc_series = - down_waves + up_wave
+    else:
+        acc_series = down_waves + up_wave
+    velocity = scipy.integrate.cumtrapz(acc_series, dx=asig.dt, initial=0, axis=1)
+    e = 0.5 * velocity * np.abs(velocity)
+    return trim_to_length(e, asig.npts, travel_times, asig.dt, trim=trim, start=start, s2s_travel_time=stt)
+
+
+def calc_cum_abs_surface_energy(asig, time_shifts, nodal=True, up_red=1, down_red=1, stt=0.0, trim=False, start=False):
+    energy = calc_surface_energy(asig, time_shifts, nodal=nodal, up_red=up_red, down_red=down_red, stt=stt,
+                                 trim=trim, start=start)
+    diff = np.zeros_like(energy)
+    diff[:, 1:] = np.diff(energy, axis=1)
+    return np.cumsum(np.abs(diff), axis=1)
+
+
+if __name__ == '__main__':
+    a = np.linspace(0, 10, 100)
+    b = np.sin(a)
+    import eqsig
+    accsig = eqsig.AccSignal(b, 0.1)
+    tshifts = np.array([0.1, 0.2])
+    ke = calc_cum_abs_surface_energy(accsig, tshifts, stt=0.3, nodal=True, up_red=1, down_red=1, trim=True)
+
diff --git a/tests/test_measures.py b/tests/test_measures.py
index 1e073c8..b3982e4 100644
--- a/tests/test_measures.py
+++ b/tests/test_measures.py
@@ -52,5 +52,11 @@ def test_calc_cyclic_amp_gm_arrays_w_power_law():
     assert np.isclose(csr_n15_series[-1], 0.96290891), (csr_n15_series[-1], 0.96290891)  # v=1.1.2
 
 
+def test_calc_unit_kinetic_energy():
+    asig = conftest.t_asig()
+    uke = im.calc_unit_kinetic_energy(asig)[-1]
+    assert np.isclose(uke, 0.4504992), uke  # version 1.1.1
+
+
 if __name__ == '__main__':
-    test_n_cyc_and_cyc_amplitude()
\ No newline at end of file
+    test_calc_unit_kinetic_energy()
diff --git a/tests/test_surface.py b/tests/test_surface.py
new file mode 100644
index 0000000..3a486a6
--- /dev/null
+++ b/tests/test_surface.py
@@ -0,0 +1,81 @@
+import numpy as np
+import eqsig
+
+from tests import conftest
+
+
+def _time_shifted_response_series(asig, times, up_red=1, down_red=1, add=False, trim=False):
+    """
+    Computes the time shifted response with reduction factors for the upward and downward motion
+
+    :param asig:
+    :param times:
+    :param up_red:
+    :param down_red:
+    :return:
+    """
+    dt = asig.dt
+    points = times / dt
+    max_points = np.max(points)
+    if trim:
+        end = asig.npts
+    else:
+        end = None
+    if hasattr(times, "__len__"):
+        if not hasattr(up_red, "__len__"):
+            up_red = up_red * np.ones_like(times)
+        if not hasattr(down_red, "__len__"):
+            down_red = down_red * np.ones_like(times)
+        energy = []
+        for tt, time in enumerate(times):
+            ss = int(points[tt])
+            extra = int(max_points - ss)
+            up_acc = list(asig.values) + [0] * ss + [0] * extra
+            down_acc = [0] * ss + list(asig.values) + [0] * extra
+            if add:
+                new_sig = np.array(up_acc) * up_red[tt] + np.array(down_acc) * down_red[tt]
+            else:
+                new_sig = np.array(up_acc) * up_red[tt] - np.array(down_acc) * down_red[tt]
+            nsig = eqsig.AccSignal(new_sig, dt)
+            ke = eqsig.im.calc_unit_kinetic_energy(nsig)
+            energy.append(np.asarray(ke[:end]))
+        return np.asarray(energy)
+    else:
+        ss = int(points)
+        extra = int(max_points - ss)
+        up_acc = list(asig.values) + [0] * ss + [0] * extra
+        down_acc = [0] * ss + list(asig.values) + [0] * extra
+        new_sig = np.array(up_acc) * up_red - np.array(down_acc) * down_red
+        nsig = eqsig.AccSignal(new_sig, dt)
+        ke = eqsig.im.calc_unit_kinetic_energy(nsig)[:end]
+        return ke
+
+
+def test_calc_cum_abs_surface_energy_w_sine_wave():
+
+    b = np.sin(np.linspace(0, 10, 100))
+    accsig = eqsig.AccSignal(b, 0.1)
+    tshifts = np.array([0.1, 0.2])
+    cases = eqsig.surface.calc_cum_abs_surface_energy(accsig, tshifts, nodal=True, up_red=1, down_red=1, stt=0.0,
+                                                      trim=True, start=False)
+    expected_cases = _time_shifted_response_series(accsig, 2 * tshifts, up_red=1, down_red=1, add=False, trim=True)
+    assert len(cases[0]) == len(b)
+    assert np.isclose(cases[0][-1], expected_cases[0][-1])
+    assert np.isclose(cases[1][-1], expected_cases[1][-1])
+
+
+def test_calc_cum_abs_surface_energy():
+    accsig = conftest.t_asig()
+    tshifts = np.array([0.01, 0.2, 0.5])
+    cases = eqsig.surface.calc_cum_abs_surface_energy(accsig, tshifts, nodal=True, up_red=1, down_red=1, stt=0.0,
+                                                      trim=True, start=False)
+    expected_cases = _time_shifted_response_series(accsig, 2 * tshifts, up_red=1, down_red=1, add=False, trim=True)
+    assert len(cases[0]) == accsig.npts
+    assert np.isclose(cases[0][-1], expected_cases[0][-1])
+    assert np.isclose(cases[1][-1], expected_cases[1][-1])
+    assert np.isclose(cases[2][-1], expected_cases[2][-1])
+
+
+if __name__ == '__main__':
+    test_calc_cum_abs_surface_energy()
+