Subtracting datasets in xarray 2024.6.0 leads to inconsistent chunks

What is your issue?

When I call groupby() on a dataset and try to subtract another dataset from the result, I get an error that says
ValueError: Object has inconsistent chunks along dimension lead. This can be fixed by calling unify_chunks().
Adding a call to unify chunks beforehand resolves the issue, but for some strange reason this chunking issue only occurs with more recent versions of xarray. When I run the same code below with xarray 2022.3.0, I can run the same code without calling unify chunks. Does anyone know what may have caused the discrepency?

Here's the relevant section of code I was running when I encountered the problem. in the snippet below, the members variable is a list of paths to netcdf files that contain the output from an ensemble of ocean models. I think the error should be reproducible with any group of netcdf files and similar operations:

ds = xarray.open_mfdataset(members, combine='nested', concat_dim='member').sortby('init') 
ensmean = ds.mean('member')
climo = ensmean.sel(init=slice('1993-01-01', '1993-12-31')).groupby('init.month').mean('init').load()
anom = model_ds.groupby('init.month') - climo

And here's the output from xarray.show_version():

INSTALLED VERSIONS
------------------
commit: None
python: 3.12.4 | packaged by conda-forge | (main, Jun 17 2024, 10:23:07) [GCC 12.3.0]
python-bits: 64
OS: Linux
OS-release: 3.10.0-1160.102.1.el7.x86_64
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: None
LANG: en_US
LOCALE: ('en_US', 'ISO8859-1')
libhdf5: 1.14.3
libnetcdf: 4.9.2

xarray: 2024.6.0
pandas: 2.2.2
numpy: 2.0.0
scipy: None
netCDF4: 1.7.1
pydap: None
h5netcdf: None
h5py: None
zarr: None
cftime: 1.6.4
nc_time_axis: None
iris: None
bottleneck: None
dask: 2024.7.1
distributed: 2024.7.1
matplotlib: None
cartopy: None
seaborn: None
numbagg: None
fsspec: 2024.6.1
cupy: None
pint: None
sparse: None
flox: None
numpy_groupies: None
setuptools: 71.0.4
pip: 24.0
conda: None
pytest: None
mypy: None
IPython: None
sphinx: None

Start by reproducing the reported open_mfdataset(), groupby(), and subtraction sequence with xarray 2024.6.0 and dask 2024.7.1, then compare it with xarray 2022.3.0. Trace where chunking changes before the inconsistent-chunks error; done means explaining the version-dependent behavior and covering it with a regression test or documented fix.