to_zarr: region not recognised as dataset dimensions

What happened:
I am trying to write into a zarr a dataset by reading each timesteps and write this time step into a zarr store. For that I prepared the dataset and filled the variable values with 0s and wrote the store with mode='a'. I then perform the operations I need on each variables and try to write the time step of the dataset.

ds.isel(time=t).to_zarr(outfile, region={"time": t})

However I received this error message:
ValueError: all keys in ``region`` are not in Dataset dimensions, got ['time'] and ['cell', 'face', 'layer', 'max_cell_node', 'max_face_nodes', 'node', 'siglay']
But

In: ds.dims
Out: Frozen(SortedKeysDict({'time': 1465, 'node': 112015, 'layer': 6, 'face': 198364, 'max_face_nodes': 3, 'cell': 991820, 'max_cell_node': 6, 'siglay': 6}))

Checking in the API, it comes from

~/.local/lib/python3.8/site-packages/xarray/backends/api.py in _validate_append_dim_and_encoding(ds_to_append, store, append_dim, region, encoding, **open_kwargs)
   1330         for k, v in region.items():
   1331             if k not in ds_to_append.dims:
-> 1332                 raise ValueError(
   1333                     f"all keys in ``region`` are not in Dataset dimensions, got "
   1334                     f"{list(region)} and {list(ds_to_append.dims)}"

What you expected to happen:
Incremently append data to the zarr store

Minimal Complete Verifiable Example:

import xarray as xr
from datetime import datetime,timedelta
import numpy as np
dt= datetime.now()
times= np.arange(dt,dt+timedelta(days=6), timedelta(hours=1))
nodesx,nodesy,layers=np.arange(10,50), np.arange(10,50)+15, np.arange(10)
ds=xr.Dataset()
ds.coords['time']=('time', times)
ds.coords['node_x']=('node', nodesx)
ds.coords['node_y']=('node', nodesy)
ds.coords['layer']=('layer', layers)
outfile='my_zarr'
varnames=['potato','banana', 'apple']
for var in varnames:
    ds[var]=(('time', 'layer', 'node'), np.zeros((len(times), len(layers),len(nodesx))))
ds.to_zarr(outfile, mode='a')
for t in range(len(times)):
    for var in varnames:
    	ds[var].isel(time=t).values += np.random.random((len(layers),len(nodesx)))
    ds.isel(time=t).to_zarr(outfile, region={"time": t})

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-25-8baf03aa01a3> in <module>
      2     for var in varnames:
      3         ds[var].isel(time=t).values += np.random.random((len(layers),len(nodesx)))
----> 4     ds.isel(time=t).to_zarr(outfile, region={"time": t})
      5 

~/.local/lib/python3.8/site-packages/xarray/core/dataset.py in to_zarr(self, store, chunk_store, mode, synchronizer, group, encoding, compute, consolidated, append_dim, region)
   1743             encoding = {}
   1744 
-> 1745         return to_zarr(
   1746             self,
   1747             store=store,

~/.local/lib/python3.8/site-packages/xarray/backends/api.py in to_zarr(dataset, store, chunk_store, mode, synchronizer, group, encoding, compute, consolidated, append_dim, region)
   1457     if mode == "a":
   1458         _validate_datatypes_for_zarr_append(dataset)
-> 1459         _validate_append_dim_and_encoding(
   1460             dataset,
   1461             store,

~/.local/lib/python3.8/site-packages/xarray/backends/api.py in _validate_append_dim_and_encoding(ds_to_append, store, append_dim, region, encoding, **open_kwargs)
   1330         for k, v in region.items():
   1331             if k not in ds_to_append.dims:
-> 1332                 raise ValueError(
   1333                     f"all keys in ``region`` are not in Dataset dimensions, got "
   1334                     f"{list(region)} and {list(ds_to_append.dims)}"

ValueError: all keys in ``region`` are not in Dataset dimensions, got ['time'] and ['layer', 'node']

Anything else we need to know?:

Environment:

Start at xarray/backends/api.py, especially _validate_append_dim_and_encoding, and reproduce the provided to_zarr loop with the stated xarray and zarr versions. Trace how Dataset.isel(time=t) changes dimensions, then add coverage for the intended region behavior while preserving validation of invalid region keys.