South facing array getting more irradiance than north facing array, even though the location is in Northern hemisphere

I am using Python and PVLib 0.10.1.

I am looking at the below example code: https://pvlib-python.readthedocs.io/en/stable/gallery/bifacial/plot_pvfactors_fixed_tilt.html#sphx-glr-gallery-bifacial-plot-pvfactors-fixed-tilt-py

"""
Fixed-Tilt Simulation with pvfactors
====================================

Modeling the irradiance on the rear side of a fixed-tilt array.
"""

# %%
# Because pvfactors was originally designed for modeling single-axis
# tracking systems, it's not necessarily obvious how to use it to model
# fixed-tilt systems correctly.
# This example shows how to model rear-side irradiance on a fixed-tilt
# array using :py:func:`pvlib.bifacial.pvfactors.pvfactors_timeseries`.
#
# .. attention::
#    To run this example, the ``solarfactors`` package (an implementation
#    of the pvfactors model) must be installed.  It can be installed with
#    either ``pip install solarfactors`` or ``pip install pvlib[optional]``,
#    which installs all of pvlib's optional dependencies.

import pandas as pd
from pvlib import location
from pvlib.bifacial.pvfactors import pvfactors_timeseries
import matplotlib.pyplot as plt
import warnings

# supressing shapely warnings that occur on import of pvfactors
warnings.filterwarnings(action='ignore', module='pvfactors')

# %%
# First, generate the usual modeling inputs:

times = pd.date_range('2021-06-21', '2021-06-22', freq='1T', tz='Etc/GMT+5')
loc = location.Location(latitude=40, longitude=-80, tz=times.tz)
sp = loc.get_solarposition(times)
cs = loc.get_clearsky(times)

# example array geometry
pvrow_height = 1
pvrow_width = 4
pitch = 10
gcr = pvrow_width / pitch
axis_azimuth = 180
albedo = 0.2

# %%
# Now the trick: since pvfactors only wants to model single-axis tracking
# arrays, we have to pretend our fixed tilt array is a single-axis tracking
# array that never rotates.  In that case, the "axis of rotation" is
# along the length of the row, with ``axis_azimuth`` 90 degrees offset from the
# fixed ``surface_azimuth``.

irrad = pvfactors_timeseries(
    solar_azimuth=sp['azimuth'],
    solar_zenith=sp['apparent_zenith'],
    surface_azimuth=180,  # south-facing array
    surface_tilt=20,
    axis_azimuth=90,  # 90 degrees off from surface_azimuth.  270 is ok too
    timestamps=times,
    dni=cs['dni'],
    dhi=cs['dhi'],
    gcr=gcr,
    pvrow_height=pvrow_height,
    pvrow_width=pvrow_width,
    albedo=albedo,
    n_pvrows=3,
    index_observed_pvrow=1
)

# turn into pandas DataFrame
irrad = pd.concat(irrad, axis=1)

irrad[['total_inc_back', 'total_abs_back']].plot()
plt.ylabel('Irradiance [W m$^{-2}$]')
plt.show()

In the example provided the default surface azimuth is to the south, and we get about 40 W / m^2 at the peak. However, when I change the surface_azimuth value to 0, i.e. facing to the North, I am getting an even higher value at around 60 W / m^2. My understanding was that I should be getting something much closer to zero as in the northern hemisphere sun will not be shining from the north.

Is my understanding of this correct, or am I doing something wrong? Why am I getting a higher value for north facing radiation than south facing radiation?

South plot:

North plot:

It is maybe worth emphasizing that you are simulating the irradiance on the rear-side (or underside) of the PV module. The typical "face towards the equator for best irradiance" rule of thumb is for front-side irradiance, since that is what aligns the module surface with the sun position the best.

Although the rear side of a PV module (typically) does not receive direct irradiance from the sun itself, it still receives diffuse irradiance from the visible portions of the ground and sky. To investigate this particular case, we have to use pvfactors itself (pip install solarfactors) instead of pvlib's convenience wrapper around pvfactors.

By extracting the irradiance components we can see that it is the ground-reflected component that is the main difference. This is because pvfactors takes ground shadows into account, and by pointing the module's front surface towards the pole, the sun can better illuminate the ground underneath the PV module. See how there is more yellow (illuminated ground) underneath the module in the azimuth=0 case.

I bet if you set the albedo to zero (effectively disabling the ground-reflected component) in your simulation, the two orientations will have negligible rear-side irradiance difference.

Here is the code to produce the above figures:

import pandas as pd
from pvlib import location
import matplotlib.pyplot as plt
import warnings
import numpy as np
from pvfactors.geometry.pvarray import OrderedPVArray
from pvfactors.engine import PVEngine

# supressing shapely warnings that occur on import of pvfactors
warnings.filterwarnings(action='ignore', module='pvfactors')

# %%

times = pd.date_range('2021-06-21', '2021-06-22', freq='1T', tz='Etc/GMT+5')
loc = location.Location(latitude=40, longitude=-80, tz=times.tz)
sp = loc.get_solarposition(times)
cs = loc.get_clearsky(times)

pvrow_height = 1
pvrow_width = 4
pitch = 10
gcr = pvrow_width / pitch
albedo = 0.2

solar_azimuth = sp['azimuth']
solar_zenith = sp['apparent_zenith']
surface_azimuth = 180
surface_tilt = 20
axis_azimuth = 90
timestamps = times
dni = cs['dni']
dhi = cs['dhi']
n_pvrows = 3
index_observed_pvrow = 1

# %%

solar_azimuth = np.array(solar_azimuth)
solar_zenith = np.array(solar_zenith)
dni = np.array(dni)
dhi = np.array(dhi)
# GH 1127, GH 1332
surface_tilt = np.full_like(solar_zenith, surface_tilt)
surface_azimuth = np.full_like(solar_zenith, surface_azimuth)

pvarray_parameters = {
    'n_pvrows': n_pvrows,
    'axis_azimuth': axis_azimuth,
    'pvrow_height': pvrow_height,
    'pvrow_width': pvrow_width,
    'gcr': gcr
}

# %%

pvarray = OrderedPVArray.init_from_dict(pvarray_parameters)
eng = PVEngine(pvarray)
eng.fit(timestamps, dni, dhi, solar_zenith, solar_azimuth, surface_tilt,
        surface_azimuth, albedo)

eng.run_full_mode()

row = pvarray.ts_pvrows[1]

data = {}
for param in row.back.all_ts_surfaces[0].params.keys():
    values = row.back.get_param_weighted(param)
    data[param] = values

df = pd.DataFrame(data, index=times)

# %%

fig, axes = plt.subplots(2, 1)
df[['direct', 'circumsolar', 'horizon', 'isotropic', 'reflection']].plot(ax=axes[0])
pvarray.plot_at_idx(len(times)//2, axes[1])
axes[0].set_title(f'{surface_azimuth=}')