diff --git a/py_wake/site/xrsite.py b/py_wake/site/xrsite.py
index ed9e56a2b455d2599af9f99a3b5f24c757d946a8..d31cc96c18e637b2f4c939702dc19af1b9717cf0 100644
--- a/py_wake/site/xrsite.py
+++ b/py_wake/site/xrsite.py
@@ -233,7 +233,9 @@ class XRSite(Site):
         lw.set_W(WS, WD, TI, ws_bins, self.use_WS_bins)
         lw.set_data_array(TI_std, 'TI_std', 'Standard deviation of turbulence intensity')
 
-        if 'time' not in lw:
+        if 'time' in lw:
+            lw['P'] = 1 / len(lw.time)
+        else:
             if 'P' in self.ds:
                 if ('ws' in self.ds.P.dims and 'ws' in lw.coords):
                     d_ws = self.ds.P.ws.values
diff --git a/py_wake/tests/test_wind_farm_models/test_enginering_wind_farm_model.py b/py_wake/tests/test_wind_farm_models/test_enginering_wind_farm_model.py
index 08f0ade0ceadfa7ddd60d802674fb07e91cab397..4a7338c938cbce08a74c6b2dbb25f5c5b0659a89 100644
--- a/py_wake/tests/test_wind_farm_models/test_enginering_wind_farm_model.py
+++ b/py_wake/tests/test_wind_farm_models/test_enginering_wind_farm_model.py
@@ -409,6 +409,18 @@ def test_time_series_override_ti():
     npt.assert_array_equal(sim_res.TI, ti)
 
 
+def test_time_series_aep():
+
+    d = np.load(os.path.dirname(examples.__file__) + "/data/time_series.npz")
+    wd, ws = [d[k][::100] for k in ['wd', 'ws']]
+    wt = V80()
+    site = Hornsrev1Site()
+    x, y = site.initial_position.T
+    wfm = NOJ(site, wt)
+    sim_res = wfm(x, y, ws=ws, wd=wd, time=True, verbose=False)
+    npt.assert_allclose(sim_res.aep().sum(), 545, atol=1)
+
+
 def test_time_series_operating():
     from py_wake.wind_turbines.power_ct_functions import PowerCtFunctionList, PowerCtTabular
     d = np.load(os.path.dirname(examples.__file__) + "/data/time_series.npz")
diff --git a/py_wake/wind_farm_models/wind_farm_model.py b/py_wake/wind_farm_models/wind_farm_model.py
index e3b0db92c34e2d72ab248cbf55620461fc9be459..42391787e061aaaf5082e26a5e2e07106d341093 100644
--- a/py_wake/wind_farm_models/wind_farm_model.py
+++ b/py_wake/wind_farm_models/wind_farm_model.py
@@ -237,6 +237,8 @@ class SimulationResult(xr.Dataset):
             return xr.DataArray(aep, [('wt', self.wt), ('wd', self.wd), ('ws', ws)])
         else:
             weighted_power = power_ilk * self.P.ilk() / norm
+        if 'time' in self and weighted_power.shape[2] == 1:
+            weighted_power = weighted_power[:, :, 0]
 
         return xr.DataArray(weighted_power * hours_pr_year * 1e-9,
                             self.Power.coords,