fixes(Mf6Splitter): multiple bug fixes and added support (#2418)

* update Mf6Splitter: Add support for groundwater energy model
* fix(split_multi_model): zero pad `model_labels` formatting string for creating exchanges
* fixes issue where "multi-models" cannot be split if there are 10 or more new model domains in the splitting mask

flopy/mf6/utils/model_splitter.py

@@ -88,6 +88,13 @@
     "lkt": "ifno",
     "mwt": "mawno",
     "uzt": "ifno",
+    # energy
+    "gwe": "cellid",
+    "ctp": "cellid",
+    "sfe": "ifno",
+    "uze": "ifno",
+
+
 }
 
 OBS_ID2_LUT = {
@@ -114,6 +121,10 @@
     },
     "mwt": None,
     "uzt": "ifno",
+    "gwe": "cellid",
+    "sfe": "ifno",
+    "uze": "ifno"
+
 }
 
 
@@ -383,13 +394,19 @@ def optimize_splitting_mask(self, nparts):
                     modflow.ModflowGwfuzf,
                     modflow.ModflowGwflak,
                     modflow.ModflowGwfhfb,
+                    modflow.ModflowGwtsft,
+                    modflow.ModflowGwtuzt,
+                    modflow.ModflowGwtlkt,
+                    modflow.ModflowGwesfe,
+                    modflow.ModflowGweuze,
+                    modflow.ModflowGwelke
                 ),
             ):
                 if isinstance(package, modflow.ModflowGwfhfb):
                     hfbs.append(package)
                     continue
 
-                if isinstance(package, modflow.ModflowGwflak):
+                if isinstance(package, (modflow.ModflowGwflak, modflow.ModflowGwtlkt, modflow.ModflowGwelke)):
                     cellids = package.connectiondata.array.cellid
                 else:
                     cellids = package.packagedata.array.cellid
@@ -401,7 +418,7 @@ def optimize_splitting_mask(self, nparts):
                 else:
                     nodes = [i[0] for i in cellids]
 
-                if isinstance(package, modflow.ModflowGwflak):
+                if isinstance(package, (modflow.ModflowGwflak, modflow.ModflowGwtlkt, modflow.ModflowGwelke)):
                     lakenos = package.connectiondata.array.ifno + 1
                     lak_array[nodes] = lakenos
                     laks += [i for i in np.unique(lakenos)]
@@ -1025,7 +1042,8 @@ def _remap_filerecords(self, item, value, mapped_data):
             "stage_filerecord",
             "obs_filerecord",
             "concentration_filerecord",
-            "ts_filerecord"
+            "ts_filerecord",
+            "temperature_filerecord"
         ):
             value = value.array
             if value is None:
@@ -1481,7 +1499,7 @@ def _remap_mvr(self, package, mapped_data):
             dict
         """
         # self._mvr_remaps = {}
-        if isinstance(package, modflow.ModflowGwtmvt):
+        if isinstance(package, (modflow.ModflowGwtmvt, modflow.ModflowGwemve)):
             return mapped_data
 
         perioddata = package.perioddata.data
@@ -2853,6 +2871,8 @@ def _remap_package(self, package, ismvr=False):
                 modflow.ModflowGwfdisu,
                 modflow.ModflowGwtdis,
                 modflow.ModflowGwtdisu,
+                modflow.ModflowGwedis,
+                modflow.ModflowGwedisu
             ),
         ):
             for item, value in package.__dict__.items():
@@ -2921,31 +2941,31 @@ def _remap_package(self, package, ismvr=False):
             mapped_data = self._remap_buy(package, mapped_data)
 
         elif isinstance(
-            package, (modflow.ModflowGwfuzf, modflow.ModflowGwtuzt)
+            package, (modflow.ModflowGwfuzf, modflow.ModflowGwtuzt, modflow.ModflowGweuze)
         ):
             mapped_data = self._remap_uzf(package, mapped_data)
 
         elif isinstance(
-            package, (modflow.ModflowGwfmaw, modflow.ModflowGwtmwt)
+            package, (modflow.ModflowGwfmaw, modflow.ModflowGwtmwt, modflow.ModflowGwemwe)
         ):
             mapped_data = self._remap_maw(package, mapped_data)
 
         elif ismvr:
             self._remap_mvr(package, mapped_data)
 
         elif isinstance(
-            package, (modflow.ModflowGwfmvr, modflow.ModflowGwtmvt)
+            package, (modflow.ModflowGwfmvr, modflow.ModflowGwtmvt, modflow.ModflowGwemve)
         ):
             self._mover = True
             return {}
 
         elif isinstance(
-            package, (modflow.ModflowGwflak, modflow.ModflowGwtlkt)
+            package, (modflow.ModflowGwflak, modflow.ModflowGwtlkt, modflow.ModflowGwelke)
         ):
             mapped_data = self._remap_lak(package, mapped_data)
 
         elif isinstance(
-            package, (modflow.ModflowGwfsfr, modflow.ModflowGwtsft)
+            package, (modflow.ModflowGwfsfr, modflow.ModflowGwtsft, modflow.ModflowGwesfe)
         ):
             mapped_data = self._remap_sfr(package, mapped_data)
 
@@ -3390,12 +3410,14 @@ def create_multi_model_exchanges(self, mname0, mname1):
 
         Parameters
         ----------
-        mname0 :
-        mname1 :
+        mname0 : str
+            model name of first model in exchange
+        mname1 : str
+            model name of second model in exchange
 
         Returns
         -------
-
+            None
         """
         exchange_classes = {
             "gwfgwt": modflow.ModflowGwfgwt,
@@ -3555,7 +3577,7 @@ def split_multi_model(self, array):
 
         gwf_base = model_names[0]
         model_labels = [
-            f"{i :{self._fdigits}d}" for i in sorted(np.unique(array))
+            f"{i :0{self._fdigits}d}" for i in sorted(np.unique(array))
         ]
 
         self._multimodel_exchange_gwf_names = {