I think you should be able to provide either a 2D array of shape (nlay, ncpl) or a 3d array of shape (nlay, 1, ncpl). Do you see an error with either of these?

The docstring does make it seem like it must be a structured grid cell index.

I provide 2d array (nlay, ncpl) and it gives error saying 3d array required. I will try doing (nlay, 1, ncpl) and can give an update if it works

Update: I have tried both. it does not work and gives the same error. I am attaching the code i am runing and the error. When i remove the stopzones, zones and zonedata option in the code the model run just fine.

npf = gwf.get_package("npf")
##get hydraulic conductivity array
k_horizontal = npf.k.array

##Store vertical hydraulic conductivity (K33) as a variable
k_vertical = npf.k33.array
#defining zones starting for my 12 layers [100,101,102,103,104,105,106,107,108,109,110 111]
layer_vals = np.arange(100, 100 + nlay).reshape(nlay, 1)

##broadcast to full array
zones = np.broadcast_to(layer_vals, (nlay, ncpl)) #resultant array (nlay,ncpl)
zones_ = np.expand_dims(zones, axis=1) #resultant array (nlay,1, ncpl)
aquifer_stopzones=[100,101,103,105,107,109,111]

pdata = ParticleData(partlocs=cellid_loc_iol,
particleids=particl_id_iol, structured=False,
localx=x_loc_iol, localy=y_loc_iol, localz=z_loc_iol, drape=1)
pg = ParticleGroup(particlegroupname="my_particles", particledata=pdata, filename="sgo.sloc")

mp7_path=r"C:\als\ALS_092125\mp7.exe"
mp = flopy.modpath.Modpath7(
modelname='mp',
flowmodel=gwf,
exe_name=mp7_path,
model_ws=sim_ws)
mpbas = flopy.modpath.Modpath7Bas(mp, porosity=0.3)
for ij in aquifer_stopzones:
mpsim = flopy.modpath.Modpath7Sim(
mp,endpointfilename='zone'+str(ij)+'.mpend',
simulationtype="endpoint",
trackingdirection="forward",
weaksinkoption="pass_through",
weaksourceoption="pass_through",
budgetoutputoption="summary",
zones=zones_,
referencetime=[0, 0, 0.0],
stoptimeoption="extend",
zonedataoption ='on',
stopzone=ij,
particlegroups=pg)
#write modpath datasets
mp.write_input()

ERROR*********
Traceback (most recent call last):

File "", line 41, in
particlegroups=pg)

File "C:\anaconda\envs\ropen\lib\site-packages\flopy\modpath\mp7sim.py", line 531, in init
locat=self.unit_number[0],

File "C:\anaconda\envs\ropen\lib\site-packages\flopy\utils\util_array.py", line 560, in init
"Util3d: expected 3 dimensions, found shape {0}".format(shape)

ValueError: Util3d: expected 3 dimensions, found shape (12, 218793)

zones_.shape
Out[15]: (12, 1, 218793)

It looks like you have a pretty old version of flopy. This issue seems to have been fixed in #1415 a few years ago. Does updating solve it?