diff --git a/Bugs_and_Todo.md b/Bugs_and_Todo.md
new file mode 100644
index 0000000..e030978
--- /dev/null
+++ b/Bugs_and_Todo.md
@@ -0,0 +1,9 @@
+## Bugs
+- [ ] HydroModel cannot used as a component in another model
+
+## TODO
+
+- [ ] Add more examples
+- [ ] Using the NamedTupleAdaptor, and remove the convert_to_ntp kwargs in component running
+- [ ] Write more error messages
+- [ ] macro building
\ No newline at end of file
diff --git a/Manifest.toml b/Manifest.toml
index dd5f114..6c399a9 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.10.4"
 manifest_format = "2.0"
-project_hash = "eec1c26075e9afae150ae4fb562f8f2fdea39275"
+project_hash = "de7e1a3bf89c00139c5d6afa1937bd039d9ad985"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "72af59f5b8f09faee36b4ec48e014a79210f2f4f"
@@ -219,6 +219,12 @@ git-tree-sha1 = "0157e592151e39fa570645e2b2debcdfb8a0f112"
 uuid = "00ebfdb7-1f24-5e51-bd34-a7502290713f"
 version = "3.4.3"
 
+[[deps.CSV]]
+deps = ["CodecZlib", "Dates", "FilePathsBase", "InlineStrings", "Mmap", "Parsers", "PooledArrays", "PrecompileTools", "SentinelArrays", "Tables", "Unicode", "WeakRefStrings", "WorkerUtilities"]
+git-tree-sha1 = "deddd8725e5e1cc49ee205a1964256043720a6c3"
+uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+version = "0.10.15"
+
 [[deps.Cassette]]
 git-tree-sha1 = "f8764df8d9d2aec2812f009a1ac39e46c33354b8"
 uuid = "7057c7e9-c182-5462-911a-8362d720325c"
@@ -246,6 +252,12 @@ git-tree-sha1 = "05ba0d07cd4fd8b7a39541e31a7b0254704ea581"
 uuid = "fb6a15b2-703c-40df-9091-08a04967cfa9"
 version = "0.1.13"
 
+[[deps.CodecZlib]]
+deps = ["TranscodingStreams", "Zlib_jll"]
+git-tree-sha1 = "bce6804e5e6044c6daab27bb533d1295e4a2e759"
+uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
+version = "0.7.6"
+
 [[deps.Combinatorics]]
 git-tree-sha1 = "08c8b6831dc00bfea825826be0bc8336fc369860"
 uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
@@ -684,6 +696,17 @@ version = "1.1.1"
     ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
     Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
 
+[[deps.FilePathsBase]]
+deps = ["Compat", "Dates"]
+git-tree-sha1 = "7878ff7172a8e6beedd1dea14bd27c3c6340d361"
+uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
+version = "0.9.22"
+weakdeps = ["Mmap", "Test"]
+
+    [deps.FilePathsBase.extensions]
+    FilePathsBaseMmapExt = "Mmap"
+    FilePathsBaseTestExt = "Test"
+
 [[deps.FileWatching]]
 uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
 
@@ -2523,6 +2546,11 @@ weakdeps = ["PDMats"]
     [deps.Tracker.extensions]
     TrackerPDMatsExt = "PDMats"
 
+[[deps.TranscodingStreams]]
+git-tree-sha1 = "0c45878dcfdcfa8480052b6ab162cdd138781742"
+uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
+version = "0.11.3"
+
 [[deps.TriangularSolve]]
 deps = ["CloseOpenIntervals", "IfElse", "LayoutPointers", "LinearAlgebra", "LoopVectorization", "Polyester", "Static", "VectorizationBase"]
 git-tree-sha1 = "be986ad9dac14888ba338c2554dcfec6939e1393"
@@ -2597,6 +2625,12 @@ git-tree-sha1 = "8351f8d73d7e880bfc042a8b6922684ebeafb35c"
 uuid = "19fa3120-7c27-5ec5-8db8-b0b0aa330d6f"
 version = "0.2.0"
 
+[[deps.WeakRefStrings]]
+deps = ["DataAPI", "InlineStrings", "Parsers"]
+git-tree-sha1 = "b1be2855ed9ed8eac54e5caff2afcdb442d52c23"
+uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5"
+version = "1.4.2"
+
 [[deps.WeightInitializers]]
 deps = ["ArgCheck", "ConcreteStructs", "GPUArraysCore", "LinearAlgebra", "Random", "SpecialFunctions", "Statistics"]
 git-tree-sha1 = "0b935265795ccccf12bc89e693b6380934451780"
@@ -2619,6 +2653,11 @@ version = "1.0.4"
     Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
     oneAPI = "8f75cd03-7ff8-4ecb-9b8f-daf728133b1b"
 
+[[deps.WorkerUtilities]]
+git-tree-sha1 = "cd1659ba0d57b71a464a29e64dbc67cfe83d54e7"
+uuid = "76eceee3-57b5-4d4a-8e66-0e911cebbf60"
+version = "1.6.1"
+
 [[deps.Zlib_jll]]
 deps = ["Libdl"]
 uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
diff --git a/Project.toml b/Project.toml
index b94a4d7..4089b45 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,6 +5,7 @@ version = "0.1.0"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
+CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
@@ -43,9 +44,9 @@ Accessors = "0.1"
 Aqua = "0.8"
 CSV = "0.10"
 ComponentArrays = "0.15"
-Dates = "1"
 DataFrames = "1"
 DataInterpolations = "6"
+Dates = "1"
 ForwardDiff = "0.10"
 Graphs = "1"
 Integrals = "4"
@@ -67,15 +68,15 @@ Reexport = "1"
 RuntimeGeneratedFunctions = "0.5"
 SciMLBase = "2"
 SciMLSensitivity = "7"
-Statistics = "1"
 SparseArrays = "1"
 StableRNGs = "1"
+Statistics = "1"
 SymbolicUtils = "3"
 Symbolics = "6"
+TOML = "1"
+Test = "1"
 Zygote = "0.6"
 julia = "1.10"
-Test = "1"
-TOML = "1"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
diff --git a/README.md b/README.md
index 303774e..d3d1611 100644
--- a/README.md
+++ b/README.md
@@ -35,3 +35,5 @@ Pkg.add("HydroModels")
 ## Contributing
 
 HydroModels.jl is an open-source project available on Github. We welcome contributions from the community to help advance deep learning applications in hydrology.
+
+
diff --git a/docs/Manifest.toml b/docs/Manifest.toml
index 56426be..41e6ca4 100644
--- a/docs/Manifest.toml
+++ b/docs/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.10.4"
 manifest_format = "2.0"
-project_hash = "ad4b8e8aabb6dad13d819587f280889b1bd2f545"
+project_hash = "620786bb65c2992bc73165cb134818fb2070b24e"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "72af59f5b8f09faee36b4ec48e014a79210f2f4f"
@@ -1056,7 +1056,7 @@ uuid = "e33a78d0-f292-5ffc-b300-72abe9b543c8"
 version = "2.11.2+1"
 
 [[deps.HydroModels]]
-deps = ["Accessors", "ComponentArrays", "DataFrames", "DataInterpolations", "Dates", "ForwardDiff", "Graphs", "Integrals", "IterTools", "LinearAlgebra", "LinearSolve", "Lux", "LuxCore", "ModelingToolkit", "NNlib", "NamedTupleTools", "Optimization", "OptimizationBBO", "OptimizationOptimisers", "OrdinaryDiffEq", "ProgressMeter", "Random", "Reexport", "RuntimeGeneratedFunctions", "SciMLBase", "SciMLSensitivity", "SparseArrays", "StableRNGs", "Statistics", "SymbolicUtils", "Symbolics", "TOML", "Zygote"]
+deps = ["Accessors", "CSV", "ComponentArrays", "DataFrames", "DataInterpolations", "Dates", "ForwardDiff", "Graphs", "Integrals", "IterTools", "LinearAlgebra", "LinearSolve", "Lux", "LuxCore", "ModelingToolkit", "NNlib", "NamedTupleTools", "Optimization", "OptimizationBBO", "OptimizationOptimisers", "OrdinaryDiffEq", "ProgressMeter", "Random", "Reexport", "RuntimeGeneratedFunctions", "SciMLBase", "SciMLSensitivity", "SparseArrays", "StableRNGs", "Statistics", "SymbolicUtils", "Symbolics", "TOML", "Zygote"]
 path = "E:\\JlCode\\HydroModels"
 uuid = "7e3cde01-c141-467b-bff6-5350a0af19fc"
 version = "0.1.0"
diff --git a/docs/Project.toml b/docs/Project.toml
index 9e6711a..86a2689 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -12,6 +12,7 @@ Lux = "b2108857-7c20-44ae-9111-449ecde12c47"
 LuxCore = "bb33d45b-7691-41d6-9220-0943567d0623"
 Measures = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+NamedTupleTools = "d9ec5142-1e00-5aa0-9d6a-321866360f50"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationBBO = "3e6eede4-6085-4f62-9a71-46d9bc1eb92b"
 OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1"
@@ -21,3 +22,4 @@ SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
+Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
diff --git a/docs/docs.rar b/docs/docs.rar
deleted file mode 100644
index 75b35c1..0000000
Binary files a/docs/docs.rar and /dev/null differ
diff --git a/docs/make.jl b/docs/make.jl
index 727bfcb..1f7e8f0 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -28,7 +28,11 @@ makedocs(
         "Tutorials" => [
             "Run a Bucket Model" => "tutorials/run_a_bucket.md",
             "Run ExpHydro Model" => "tutorials/run_a_exphydro_model.md"
-        ]
+        ],
+        "Extent Content" => [
+            "Modeling Framework Comparisons" => "docs/src/extent/framework comparision - en.md",
+        ],
+        
     ],
     # 其他选项
     checkdocs = :none,
diff --git a/docs/src/extent/framework comparision - en.md b/docs/src/extent/framework comparision - en.md
new file mode 100644
index 0000000..b0e820c
--- /dev/null
+++ b/docs/src/extent/framework comparision - en.md	
@@ -0,0 +1,142 @@
+# Comparing Programming Styles for Custom Models in superflexpy, MARRMoT, and HydroModels.jl
+
+This content corresponds to Section 4.1 of the paper, where we discuss the differences in programming styles for custom model implementation across these three software libraries. Using the GR4J model implementation as an example, we compare the programming styles for customizing model computation formulas.
+
+## [Superflexpy](https://github.com/dalmo1991/superflexPy/tree/master)
+
+Superflexpy is a Python implementation based on the SUPERFLEX concept. For implementing hydrological model computations, this framework typically uses Python syntax to express model calculation formulas, as shown below.
+
+```python
+@staticmethod
+def _flux_function_python(S, S0, ind, P, x1, alpha, beta, ni, PET, dt):
+    if ind is None:
+        return (
+            [
+                P * (1 - (S / x1) ** alpha), # Ps
+                -PET * (2 * (S / x1) - (S / x1) ** alpha), # Evaporation
+                -((x1 ** (1 - beta)) / ((beta - 1) * dt)) * (ni ** (beta - 1)) * (S**beta), # Perc
+            ],
+            0.0,
+            S0 + P * (1 - (S / x1) ** alpha) * dt,
+        )
+    else:
+        return (
+            [
+                P[ind] * (1 - (S / x1[ind]) ** alpha[ind]), # Ps
+                -PET[ind] * (2 * (S / x1[ind]) - (S / x1[ind]) ** alpha[ind]), # Evaporation
+                -((x1[ind] ** (1 - beta[ind])) / ((beta[ind] - 1) * dt[ind]))
+                * (ni[ind] ** (beta[ind] - 1))
+                * (S ** beta[ind]), # Perc
+            ],
+            0.0,
+            S0 + P[ind] * (1 - (S / x1[ind]) ** alpha[ind]) * dt[ind],
+            [
+                -(P[ind] * alpha[ind] / x1[ind]) * ((S / x1[ind]) ** (alpha[ind] - 1)),
+                -(PET[ind] / x1[ind]) * (2 - alpha[ind] * ((S / x1[ind]) ** (alpha[ind] - 1))),
+                -beta[ind]
+                * ((x1[ind] ** (1 - beta[ind])) / ((beta[ind] - 1) * dt[ind]))
+                * (ni[ind] ** (beta[ind] - 1))
+                * (S ** (beta[ind] - 1)),
+            ],
+        )
+
+@staticmethod
+@nb.jit(
+    "Tuple((UniTuple(f8, 3), f8, f8, UniTuple(f8, 3)))"
+    "(optional(f8), f8, i4, f8[:], f8[:], f8[:], f8[:], f8[:], f8[:], f8[:])",
+    nopython=True,
+)
+def _flux_function_numba(S, S0, ind, P, x1, alpha, beta, ni, PET, dt):
+    return (
+        (
+            P[ind] * (1 - (S / x1[ind]) ** alpha[ind]),  # Ps
+            -PET[ind] * (2 * (S / x1[ind]) - (S / x1[ind]) ** alpha[ind]),  # Evaporation
+            -((x1[ind] ** (1 - beta[ind])) / ((beta[ind] - 1) * dt[ind]))
+            * (ni[ind] ** (beta[ind] - 1))
+            * (S ** beta[ind]),  # Perc
+        ),
+        0.0,
+        S0 + P[ind] * (1 - (S / x1[ind]) ** alpha[ind]) * dt[ind],
+        (
+            -(P[ind] * alpha[ind] / x1[ind]) * ((S / x1[ind]) ** (alpha[ind] - 1)),
+            -(PET[ind] / x1[ind]) * (2 - alpha[ind] * ((S / x1[ind]) ** (alpha[ind] - 1))),
+            -beta[ind]
+            * ((x1[ind] ** (1 - beta[ind])) / ((beta[ind] - 1) * dt[ind]))
+            * (ni[ind] ** (beta[ind] - 1))
+            * (S ** (beta[ind] - 1)),
+        ),
+    )
+```
+
+As shown, when implementing the GR4J production module calculations [origin code](https://github.com/dalmo1991/superflexPy/blob/master/superflexpy/implementation/elements/gr4j.py), the superflexpy framework requires separate construction of numpy and numba computation code. The function must specify hydrological fluxes `Ps`, `Evaporation`, `Perc`, and state variable balance equations. This approach couples multiple computational process implementations in a single code block, resulting in reduced code reusability (only reusable at the module level) and increased code volume and redundancy.
+
+## [MARRMoT](https://github.com/wknoben/MARRMoT)
+
+MARRMoT is a hydrological modeling framework developed in MATLAB. This framework constructs model types separately based on different hydrological models, storing the calculation formulas and state change equations involved. The simplified implementation code for the GR4J model in this framework is shown below.
+
+
+```matlab
+function [dS, fluxes] = model_fun(obj, S)
+    % definitions
+    ...
+    % fluxes functions
+    flux_pn   = max(P-Ep,0);
+    flux_en   = max(Ep-P,0);
+    flux_ef   = P - flux_pn;
+    flux_ps   = saturation_4(S1,x1,flux_pn);
+    flux_es   = evap_11(S1,x1,flux_en);
+    flux_perc = percolation_3(S1,x1);
+    flux_q9   = route(.9.*(flux_pn - flux_ps + flux_perc), uh_q9);
+    flux_q1   = route(.1.*(flux_pn - flux_ps + flux_perc), uh_q1);
+    flux_fr   = recharge_2(3.5,S2,x3,x2);
+    flux_fq   = flux_fr;
+    flux_qr   = baseflow_3(S2,x3);
+    flux_qt   = flux_qr + max(flux_q1 + flux_fq,0);
+    % this flux is not included in original MARRMoT,
+    % but it is useful to calculate the water balance
+    flux_ex = flux_fr + max(flux_q1 + flux_fq,0) - flux_q1;      
+    % stores ODEs
+    dS1 = flux_ps - flux_es - flux_perc;
+    dS2 = flux_q9 + flux_fr - flux_qr;
+    % outputs
+    dS = [dS1 dS2];
+    fluxes = [flux_pn,   flux_en, flux_ef, flux_ps, flux_es,...
+                flux_perc, flux_q9, flux_q1, flux_fr, flux_fq,...
+                flux_qr,   flux_qt, flux_ex];
+end
+```
+
+As evident, the MARRMoT framework expresses each Flux calculation formula through functions, allowing users to directly call functions for calculations (according to provided documentation). This approach significantly simplifies code writing requirements compared to superflexpy. However, MARRMoT has poor support for module reusability and cannot be reused like superflexpy. Additionally, combining variable assignments, flux calculations, and state change equations in one function results in poor decoupling of computational content.
+
+## HydroModels.jl
+
+```julia
+#* define the production store
+prod_funcs = [
+    HydroFlux([prcp, ep] => [pn, en], exprs=[prcp - min(prcp, ep), ep - min(prcp, ep)]),
+    HydroFlux([pn, soilwater] => [ps], [x1], exprs=[max(0.0, pn * (1 - (soilwater / x1)^2))]),
+    HydroFlux([en, soilwater] => [es], [x1], exprs=[en * (2 * soilwater / x1 - (soilwater / x1)^2)]),
+    HydroFlux([soilwater] => [perc], [x1], exprs=[((x1)^(-4)) / 4 * ((4 / 9)^(4)) * (soilwater^5)]),
+    HydroFlux([pn, ps, perc] => [pr], exprs=[pn - ps + perc]),
+    HydroFlux([pr] => [slowflow, fastflow], exprs=[0.9 * pr, 0.1 * pr]),
+]
+prod_dfuncs = [StateFlux([ps] => [es, perc], soilwater)]
+prod_ele = HydroBucket(name=:gr4j_prod, funcs=prod_funcs, dfuncs=prod_dfuncs)
+
+#* uh function
+uh_flux_1 = HydroModels.UnitHydrograph([slowflow] => [slowflow_routed], x4, uhfunc=HydroModels.UHFunction(:UH_1_HALF), solvetype=:SPARSE)
+uh_flux_2 = HydroModels.UnitHydrograph([fastflow] => [fastflow_routed], x4, uhfunc=HydroModels.UHFunction(:UH_2_FULL), solvetype=:SPARSE)
+
+#* define the routing store
+rst_funcs = [
+    HydroFlux([routingstore] => [exch], [x2, x3], exprs=[x2 * abs(routingstore / x3)^3.5]),
+    HydroFlux([routingstore, slowflow, exch] => [routedflow], [x3], exprs=[x3^(-4) / 4 * (routingstore + slowflow + exch)^5]),
+    HydroFlux([routedflow, fastflow_routed, exch] => [flow], exprs=[routedflow + max(fastflow_routed + exch, 0.0)]),
+]
+rst_dfuncs = [StateFlux([exch, slowflow] => [routedflow], routingstore)]
+rst_ele = HydroBucket(name=:gr4j_rst, funcs=rst_funcs, dfuncs=rst_dfuncs)
+
+gr4j_model = HydroModel(name=:gr4j, components=[prod_ele, uh_flux_1, uh_flux_2, rst_ele])
+```
+
+Returning to the GR4J model implementation in HydroModels.jl framework, we sequentially define the production store, unit hydrograph, and routing store. In this process, the model can achieve reusability while maintaining variable correspondence. The construction process includes instantiation of classes such as HydroFlux, StateFlux, UnitHydrograph, HydroModel, and HydroBucket. Unlike the previous two frameworks that use functional construction, these class instantiations express calculation formulas through symbolic programming. After model construction, calculations are performed through input data and parameters, demonstrating a design philosophy that completely decouples model structure and formulas from parameters and data, achieving high reusability.
\ No newline at end of file
diff --git a/docs/src/extent/framework comparision - zh.md b/docs/src/extent/framework comparision - zh.md
new file mode 100644
index 0000000..f588674
--- /dev/null
+++ b/docs/src/extent/framework comparision - zh.md	
@@ -0,0 +1,141 @@
+# 比较superflexpy, MARRMoT 和 HydroModels.jl三种库对于自定义模型的编程风格
+
+本内容对应了论文的4.1节的内容,我们在这里希望能够讨论三种软件库对于模型计算公式的自定义的编程风格的差异. 我们以GR4J模型的实习过程为例,对模型计算公式的自定义的编程风格进行比较.
+
+## [Superflexpy](https://github.com/dalmo1991/superflexPy/tree/master)
+
+Superflexpy是基于SUPERFLEX思想下的python代码实现, 针对于水文模型中一些计算功能的实现, 该框架一般而言是直接使用python语法表述模型的计算公式, 如下所示.
+
+```python
+@staticmethod
+def _flux_function_python(S, S0, ind, P, x1, alpha, beta, ni, PET, dt):
+    if ind is None:
+        return (
+            [
+                P * (1 - (S / x1) ** alpha), # Ps
+                -PET * (2 * (S / x1) - (S / x1) ** alpha), # Evaporation
+                -((x1 ** (1 - beta)) / ((beta - 1) * dt)) * (ni ** (beta - 1)) * (S**beta), # Perc
+            ],
+            0.0,
+            S0 + P * (1 - (S / x1) ** alpha) * dt,
+        )
+    else:
+        return (
+            [
+                P[ind] * (1 - (S / x1[ind]) ** alpha[ind]), # Ps
+                -PET[ind] * (2 * (S / x1[ind]) - (S / x1[ind]) ** alpha[ind]), # Evaporation
+                -((x1[ind] ** (1 - beta[ind])) / ((beta[ind] - 1) * dt[ind]))
+                * (ni[ind] ** (beta[ind] - 1))
+                * (S ** beta[ind]), # Perc
+            ],
+            0.0,
+            S0 + P[ind] * (1 - (S / x1[ind]) ** alpha[ind]) * dt[ind],
+            [
+                -(P[ind] * alpha[ind] / x1[ind]) * ((S / x1[ind]) ** (alpha[ind] - 1)),
+                -(PET[ind] / x1[ind]) * (2 - alpha[ind] * ((S / x1[ind]) ** (alpha[ind] - 1))),
+                -beta[ind]
+                * ((x1[ind] ** (1 - beta[ind])) / ((beta[ind] - 1) * dt[ind]))
+                * (ni[ind] ** (beta[ind] - 1))
+                * (S ** (beta[ind] - 1)),
+            ],
+        )
+
+@staticmethod
+@nb.jit(
+    "Tuple((UniTuple(f8, 3), f8, f8, UniTuple(f8, 3)))"
+    "(optional(f8), f8, i4, f8[:], f8[:], f8[:], f8[:], f8[:], f8[:], f8[:])",
+    nopython=True,
+)
+def _flux_function_numba(S, S0, ind, P, x1, alpha, beta, ni, PET, dt):
+    return (
+        (
+            P[ind] * (1 - (S / x1[ind]) ** alpha[ind]),  # Ps
+            -PET[ind] * (2 * (S / x1[ind]) - (S / x1[ind]) ** alpha[ind]),  # Evaporation
+            -((x1[ind] ** (1 - beta[ind])) / ((beta[ind] - 1) * dt[ind]))
+            * (ni[ind] ** (beta[ind] - 1))
+            * (S ** beta[ind]),  # Perc
+        ),
+        0.0,
+        S0 + P[ind] * (1 - (S / x1[ind]) ** alpha[ind]) * dt[ind],
+        (
+            -(P[ind] * alpha[ind] / x1[ind]) * ((S / x1[ind]) ** (alpha[ind] - 1)),
+            -(PET[ind] / x1[ind]) * (2 - alpha[ind] * ((S / x1[ind]) ** (alpha[ind] - 1))),
+            -beta[ind]
+            * ((x1[ind] ** (1 - beta[ind])) / ((beta[ind] - 1) * dt[ind]))
+            * (ni[ind] ** (beta[ind] - 1))
+            * (S ** (beta[ind] - 1)),
+        ),
+    )
+```
+
+可以看出, superflexpy框架在实现GR4J的production模块计算时[origin code](https://github.com/dalmo1991/superflexPy/blob/master/superflexpy/implementation/elements/gr4j.py),不仅需要分别构建numpy和numba的计算代码,在函数中需要给出模型中涉及的水文通量`Ps`,`Evaporation`,`Perc`,同时还要给出状态变量的平衡方程,一个代码中需要耦合多个计算过程的实现代码,不仅导致代码的复用的灵活性较差(仅能在模块层面进行复用),模型编写的代码量和重复度也较高.
+
+## [MARRMoT](https://github.com/wknoben/MARRMoT)
+
+MARRMoT是一个由MATLAB语言开发的水文模型框架,该框架根据不同水文模型会分别构建模型类型,储存模型涉及的计算公式和状态变化方程,GR4J模型在该框架下的实现代码(简化了变量命名)如下所示.
+
+```matlab
+function [dS, fluxes] = model_fun(obj, S)
+    % definitions
+    ...
+    % fluxes functions
+    flux_pn   = max(P-Ep,0);
+    flux_en   = max(Ep-P,0);
+    flux_ef   = P - flux_pn;
+    flux_ps   = saturation_4(S1,x1,flux_pn);
+    flux_es   = evap_11(S1,x1,flux_en);
+    flux_perc = percolation_3(S1,x1);
+    flux_q9   = route(.9.*(flux_pn - flux_ps + flux_perc), uh_q9);
+    flux_q1   = route(.1.*(flux_pn - flux_ps + flux_perc), uh_q1);
+    flux_fr   = recharge_2(3.5,S2,x3,x2);
+    flux_fq   = flux_fr;
+    flux_qr   = baseflow_3(S2,x3);
+    flux_qt   = flux_qr + max(flux_q1 + flux_fq,0);
+    % this flux is not included in original MARRMoT,
+    % but it is useful to calculate the water balance
+    flux_ex = flux_fr + max(flux_q1 + flux_fq,0) - flux_q1;      
+    % stores ODEs
+    dS1 = flux_ps - flux_es - flux_perc;
+    dS2 = flux_q9 + flux_fr - flux_qr;
+    % outputs
+    dS = [dS1 dS2];
+    fluxes = [flux_pn,   flux_en, flux_ef, flux_ps, flux_es,...
+                flux_perc, flux_q9, flux_q1, flux_fr, flux_fq,...
+                flux_qr,   flux_qt, flux_ex];
+end
+```
+
+可以看出,MARRMoT框架将每一个Flux计算公式由函数进行表达,因此使用者可以直接调取函数进行计算(可以根据提供的文档进行调用),这种方式相比superflexpy的确简化了很多代码编写需求. 然而MARRMoT对于模块的重复率支持性较差,不能够像superflexpy那样进行复用. 同时将变量赋值,通量计算和状态变化方程以一个函数包括,计算内容的解耦性较差.
+
+## HydroModels.jl
+
+```julia
+#* define the production store
+prod_funcs = [
+    HydroFlux([prcp, ep] => [pn, en], exprs=[prcp - min(prcp, ep), ep - min(prcp, ep)]),
+    HydroFlux([pn, soilwater] => [ps], [x1], exprs=[max(0.0, pn * (1 - (soilwater / x1)^2))]),
+    HydroFlux([en, soilwater] => [es], [x1], exprs=[en * (2 * soilwater / x1 - (soilwater / x1)^2)]),
+    HydroFlux([soilwater] => [perc], [x1], exprs=[((x1)^(-4)) / 4 * ((4 / 9)^(4)) * (soilwater^5)]),
+    HydroFlux([pn, ps, perc] => [pr], exprs=[pn - ps + perc]),
+    HydroFlux([pr] => [slowflow, fastflow], exprs=[0.9 * pr, 0.1 * pr]),
+]
+prod_dfuncs = [StateFlux([ps] => [es, perc], soilwater)]
+prod_ele = HydroBucket(name=:gr4j_prod, funcs=prod_funcs, dfuncs=prod_dfuncs)
+
+#* uh function
+uh_flux_1 = HydroModels.UnitHydrograph([slowflow] => [slowflow_routed], x4, uhfunc=HydroModels.UHFunction(:UH_1_HALF), solvetype=:SPARSE)
+uh_flux_2 = HydroModels.UnitHydrograph([fastflow] => [fastflow_routed], x4, uhfunc=HydroModels.UHFunction(:UH_2_FULL), solvetype=:SPARSE)
+
+#* define the routing store
+rst_funcs = [
+    HydroFlux([routingstore] => [exch], [x2, x3], exprs=[x2 * abs(routingstore / x3)^3.5]),
+    HydroFlux([routingstore, slowflow, exch] => [routedflow], [x3], exprs=[x3^(-4) / 4 * (routingstore + slowflow + exch)^5]),
+    HydroFlux([routedflow, fastflow_routed, exch] => [flow], exprs=[routedflow + max(fastflow_routed + exch, 0.0)]),
+]
+rst_dfuncs = [StateFlux([exch, slowflow] => [routedflow], routingstore)]
+rst_ele = HydroBucket(name=:gr4j_rst, funcs=rst_funcs, dfuncs=rst_dfuncs)
+
+gr4j_model = HydroModel(name=:gr4j, components=[prod_ele, uh_flux_1, uh_flux_2, rst_ele])
+```
+
+回到HydroModels.jl框架中对于GR4J模型的实现, 在该代码中我们依次对production store, unit hydrograph, routing store进行了定义, 在这个过程中,在变量对应的前提下,模型能够实现可复用的能力. 在这个搭建过程中, 包括HydroFlux, StateFlux, UnitHydrograph, HydroModel, HydroBucket等类的实例化,这些类型的实例化方式不是像前两者一样使用函数的方式构建的, 而是通过符号编程的方式对计算公式进行表达, 在模型构建完成后是通过输入数据和参数进行模型计算,可以看出这种设计思路是将模型结构与公式同参数和数据完全解耦,有着极高的可复用性.
\ No newline at end of file
diff --git a/docs/src/tutorials/build_and_run_m50.md b/docs/src/tutorials/build_and_run_m50.md
new file mode 100644
index 0000000..7a57a81
--- /dev/null
+++ b/docs/src/tutorials/build_and_run_m50.md
@@ -0,0 +1 @@
+# Embedding Neural Networks in ExpHydro Model: Implementation of [M50 and M100](https://hess.copernicus.org/articles/26/5085/2022/) Model
\ No newline at end of file
diff --git a/docs/src/tutorials/run_a_exphydro_model.md b/docs/src/tutorials/run_a_exphydro_model.md
index a26c5d3..38e299a 100644
--- a/docs/src/tutorials/run_a_exphydro_model.md
+++ b/docs/src/tutorials/run_a_exphydro_model.md
@@ -1,57 +1,69 @@
-# Run A ExpHydro Model
+# ExpHydro Model Tutorial
 
 ## Overview
 
-本文档将介绍如何使用已构建的exphdyro模型用于集总式(Single HRU)的计算和分布式(Multiple HURs)的计算
+This tutorial demonstrates how to use the ExpHydro model for hydrological calculations, including:
+- Lumped calculation (Single HRU)
+- Distributed calculation (Multiple HRUs)
 
 ## Dependencies
 
-First, lets import the required packages.
+First, import the required Julia packages:
 
 ```julia
-using CSV
-using DataFrames
-using ComponentArrays
-using BenchmarkTools
-using NamedTupleTools
-using Plots
+using CSV            # Data reading
+using DataFrames     # Data processing
+using ComponentArrays # Parameter organization
+using BenchmarkTools # Performance testing
+using NamedTupleTools # Named tuple utilities
+using Plots         # Result visualization
 ```
 
 ## Model Setup
 
-接着需要对模型的参数和初始状态进行设置
-
 ### Parameter Configuration
 
-Exphydro模型包含以下6个参数,我们需要使用ComponentVector来定义这些参数
+The ExpHydro model contains 6 key parameters defined using ComponentVector:
 
 | Parameter | Value        | Description                   |
 | --------- | ------------ | ----------------------------- |
-| `f`       | 0.01674478   | Infiltration parameter        |
-| `Smax`    | 1709.461015  | Maximum soil water storage    |
-| `Qmax`    | 18.46996175  | Maximum discharge             |
-| `Df`      | 2.674548848  | Degree-day factor             |
-| `Tmax`    | 0.175739196  | Maximum temperature threshold |
-| `Tmin`    | -2.092959084 | Minimum temperature threshold |
+| `f`       | 0.01674478  | Infiltration parameter        |
+| `Smax`    | 1709.461015 | Maximum soil water storage    |
+| `Qmax`    | 18.46996175 | Maximum discharge             |
+| `Df`      | 2.674548848 | Degree-day factor             |
+| `Tmax`    | 0.175739196 | Maximum temperature threshold |
+| `Tmin`    | -2.092959084| Minimum temperature threshold |
 
 ```julia
-params = ComponentVector(f=0.01674478, Smax=1709.461015, Qmax=18.46996175,
- Df=2.674548848, Tmax=0.175739196, Tmin=-2.092959084)
+# Define model parameters
+params = ComponentVector(
+    f = 0.01674478,    # Infiltration parameter
+    Smax = 1709.461015, # Maximum soil water storage
+    Qmax = 18.46996175, # Maximum discharge
+    Df = 2.674548848,   # Degree-day factor
+    Tmax = 0.175739196, # Maximum temperature threshold
+    Tmin = -2.092959084 # Minimum temperature threshold
+)
 ```
 
 ### Initial States
 
-然后就是定义模型的初始状态,针对exphydro模型的两个计算模块分别设置其对应的状态变量`snowpack`和`soilwater`初始值
+Set initial states for the two calculation modules:
 
 ```julia
-inistates = ComponentVector(snowpack=0.0, soilwater=1303.004248)
+# Define initial states
+inistates = ComponentVector(
+    snowpack = 0.0,           # Snow accumulation
+    soilwater = 1303.004248   # Soil water content
+)
 ```
 
 ## Data Preparation
 
-The test uses hydrometeorological data from "data/exphydro/01013500.csv":
+The model uses hydrometeorological data from "data/exphydro/01013500.csv":
 
 ```julia
+# Read input data
 file_path = "data/exphydro/01013500.csv"
 data = CSV.File(file_path)
 df = DataFrame(data)
@@ -72,41 +84,41 @@ Input variables include:
 # Configure solver
 solver = HydroModels.ManualSolver()
 
-# Run model with performance benchmarking
-result = model(input, pas, config=(solver=solver, timeidx=ts), convert_to_ntp=true)
+# Run model with benchmarking
+result = model(input, params, config=(solver=solver, timeidx=ts), convert_to_ntp=true)
 
 # Visualize results
-plot(result.flow)
-plot!(df[ts, "flow(mm)"])
+plot(result.flow, label="Simulated")
+plot!(df[ts, "flow(mm)"], label="Observed")
 ```
 
-### Performance Benchmarking For Single Node (by using BenchmarkTools.jl)
+### Single Node Benchmarking
 
 ```julia
-@btime model(input, pas, config=(solver=solver, timeidx=ts), convert_to_ntp=true);
+# Performance testing using BenchmarkTools
+@btime model(input, params, config=(solver=solver, timeidx=ts), convert_to_ntp=true);
 ```
 
-### Multi-Node Testing (Optional)
+### Multi-Node Testing
 
-The code includes commented sections for multi-node testing:
+Multi-node setup for distributed computing:
 
 ```julia
 # Setup multiple nodes
-node_num = 10
+node_num = 10  # Number of nodes
 inputs = repeat([input], node_num)
 ptypes = [Symbol(:node, i) for i in 1:node_num]
 
 # Configure parameters and states for multiple nodes
 params_multi = ComponentVector(NamedTuple{Tuple(ptypes)}(repeat([params], node_num)))
-init_states_multi = ComponentVector(NamedTuple{Tuple(ptypes)}(repeat([init_states], node_num)))
-pas_multi = ComponentVector(params=params_multi, initstates=init_states_multi)
+init_states_multi = ComponentVector(NamedTuple{Tuple(ptypes)}(repeat([inistates], node_num)))
 
 # Run multi-node simulation
-results = model(inputs, pas_multi, config=(solver=solver, timeidx=ts), convert_to_ntp=true)
+results = model(inputs, params_multi, config=(solver=solver, timeidx=ts), convert_to_ntp=true)
 ```
 
-### Performance Benchmarking For Multi-Node (by using BenchmarkTools.jl)
+### Multi-Node Benchmarking
 
 ```julia
-@btime model(inputs, pas_multi, config=(solver=solver, timeidx=ts), convert_to_ntp=true);
-```
\ No newline at end of file
+# Multi-node performance testing
+@btime model(inputs, params_multi, config=(solver=solver, timeidx=ts), convert_to_ntp=true);
\ No newline at end of file
diff --git a/docs/src/tutorials/run_fluxes.md b/docs/src/tutorials/run_fluxes.md
new file mode 100644
index 0000000..e69de29
diff --git a/src/HydroModels.jl b/src/HydroModels.jl
index 4bd023a..1f968f6 100644
--- a/src/HydroModels.jl
+++ b/src/HydroModels.jl
@@ -64,6 +64,7 @@ struct HydroEquation end
 abstract type AbstractComponent end
 abstract type AbstractHydroSolver end
 abstract type AbstractHydroOptimizer end
+abstract type AbstractIOAdapter end
 abstract type AbstractHydroWrapper <: AbstractComponent end
 abstract type AbstractNeuralWrapper <: AbstractComponent end
 
@@ -92,6 +93,8 @@ include("utils/show.jl")
 include("utils/build.jl")
 include("utils/callback.jl")
 include("utils/sort.jl")
+include("utils/io.jl")
+export NamedTupleIOAdapter
 
 include("optimizer.jl")
 export BatchOptimizer, HydroOptimizer, GradOptimizer
diff --git a/src/bucket.jl b/src/bucket.jl
index ad4988b..3b71a93 100644
--- a/src/bucket.jl
+++ b/src/bucket.jl
@@ -289,14 +289,15 @@ function (ele::HydroBucket{F,D,FF,OF,M})(
     ele_output_vec = [ele.flux_func.(eachslice(fluxes[:, :, i], dims=2), param_func(pas), nn_param_func(pas), timeidx[i]) for i in axes(fluxes, 3)]
     ele_output_arr = reduce((m1, m2) -> cat(m1, m2, dims=3), [reduce(hcat, u) for u in ele_output_vec])
     #* merge state and output, if solved_states is not nothing, then cat it at the first dim
-    final_output_arr = cat(solved_states, ele_output_arr, dims=1)
+    combine_output_arr = cat(solved_states, ele_output_arr, dims=1)
 
     #* convert to NamedTuple if convert_to_ntp is true
     convert_to_ntp = get(kwargs, :convert_to_ntp, false)
     if convert_to_ntp
-        return [NamedTuple{Tuple(vcat(get_state_names(ele), get_output_names(ele)))}(eachslice(final_output_arr[:, i, :], dims=1)) for i in axes(final_output_arr, 2)]
+        output_names_tuple = Tuple(vcat(get_state_names(ele), get_output_names(ele)))
+        return [NamedTuple{output_names_tuple}(eachslice(output_arr_, dims=1)) for output_arr_ in eachslice(combine_output_arr, dims=2)]
     else
-        return final_output_arr
+        return combine_output_arr
     end
 end
 
@@ -326,7 +327,8 @@ function (ele::HydroBucket{F,D,FF,OF,M})(
     #* convert to NamedTuple if convert_to_ntp is true
     convert_to_ntp = get(kwargs, :convert_to_ntp, false)
     if convert_to_ntp
-        return [NamedTuple{Tuple(vcat(get_state_names(ele), get_output_names(ele)))}(eachslice(final_output_arr[:, i, :], dims=1)) for i in axes(final_output_arr, 2)]
+        output_names_tuple = Tuple(vcat(get_state_names(ele), get_output_names(ele)))
+        return [NamedTuple{output_names_tuple}(eachslice(output_arr_, dims=1)) for output_arr_ in eachslice(final_output_arr, dims=2)]
     else
         return final_output_arr
     end
diff --git a/src/flux.jl b/src/flux.jl
index 5b35da9..7cafdcd 100644
--- a/src/flux.jl
+++ b/src/flux.jl
@@ -118,14 +118,7 @@ function (flux::AbstractHydroFlux)(input::AbstractArray{N,3}, pas::ComponentVect
     @assert length(timeidx) == size(input, 3) "Time index length does not match the number of time steps"
 
     #* extract params and nn params
-    node_params = [pas[:params][ptype] for ptype in ptypes]
-    # 批量检查参数
-    required_params = Set(get_param_names(flux))
-    for (ptype, node_param) in zip(ptypes, node_params)
-        missing_params = setdiff(required_params, keys(node_param))
-        @assert isempty(missing_params) "Missing parameters for $ptype: $missing_params"
-    end
-    params_vec = collect([collect([params_item[pname] for pname in get_param_names(flux)]) for params_item in node_params])
+    params_vec = collect([collect([pas[:params][ptype][pname] for pname in get_param_names(flux)]) for ptype in ptypes])
 
     #* array dims: (var_names * node_names * ts_len)
     flux_output_vec = [reduce(hcat, flux.func.(eachslice(input[:, :, i], dims=2), params_vec, timeidx[i])) for i in eachindex(timeidx)]
@@ -135,9 +128,6 @@ function (flux::AbstractHydroFlux)(input::AbstractArray{N,3}, pas::ComponentVect
     else
         flux_output_arr = reduce((m1, m2) -> cat(m1, m2, dims=3), flux_output_vec)
     end
-    # flux_output_arr = mapslices(input, dims=[1,2]) do slice
-    #     reduce(hcat, flux.func.(eachcol(slice), params_vec, timeidx))
-    # end
     flux_output_arr
 end
 
diff --git a/src/model.jl b/src/model.jl
index da2a805..a559d6b 100644
--- a/src/model.jl
+++ b/src/model.jl
@@ -119,5 +119,5 @@ function (model::HydroModel)(
         fluxes = cat(fluxes, tmp_fluxes, dims=1)
     end
     convert_to_ntp = get(kwargs, :convert_to_ntp, false)
-    return convert_to_ntp ? [NamedTuple{Tuple(model.varnames)}(eachslice(fluxes[:, i, :], dims=1)) for i in axes(fluxes, 2)] : fluxes
+    return convert_to_ntp ? [NamedTuple{Tuple(model.varnames)}(eachslice(fluxes_, dims=1)) for fluxes_ in eachslice(fluxes, dims=2)] : fluxes
 end
diff --git a/src/optimizer.jl b/src/optimizer.jl
index be14ab3..63a543d 100644
--- a/src/optimizer.jl
+++ b/src/optimizer.jl
@@ -74,7 +74,7 @@ function (opt::HydroOptimizer{C,S})(
     target::Vector;
     tunable_pas::ComponentVector,
     const_pas::ComponentVector,
-    config::Vector{<:NamedTuple}=fill(NamedTuple(), length(input)),
+    config::Vector=fill(NamedTuple(), length(input)),
     kwargs...
 ) where {C,S}
     lb = get(kwargs, :lb, zeros(length(tunable_pas)))
@@ -89,6 +89,7 @@ function (opt::HydroOptimizer{C,S})(
     prob_args = (input, target, config, run_kwargs, tunable_axes, default_model_pas)
     #* Constructing and solving optimization problems
     optf = Optimization.OptimizationFunction(opt.objective_func)
+    @info "The size of tunable parameters is $(length(tunable_pas))"
     optprob = Optimization.OptimizationProblem(optf, collect(tunable_pas), prob_args, lb=lb, ub=ub)
     sol = Optimization.solve(optprob, opt.solve_alg, callback=callback_func, maxiters=opt.maxiters)
     opt_pas = update_ca(default_model_pas, ComponentVector(sol.u, tunable_axes))
@@ -105,7 +106,7 @@ function (opt::GradOptimizer{C,S})(
     target::Vector;
     tunable_pas::ComponentVector,
     const_pas::ComponentVector,
-    config::Vector{<:NamedTuple}=fill(NamedTuple(), length(input)),
+    config::Vector=fill(NamedTuple(), length(input)),
     kwargs...
 ) where {C,S}
     run_kwargs = get(kwargs, :run_kwargs, (convert_to_ntp=true,))
@@ -117,6 +118,7 @@ function (opt::GradOptimizer{C,S})(
     prob_args = (input, target, config, run_kwargs, tunable_axes, default_model_pas)
     #* Constructing and solving optimization problems
     optf = Optimization.OptimizationFunction(opt.objective_func, opt.adtype)
+    @info "The size of tunable parameters is $(length(tunable_pas))"
     optprob = Optimization.OptimizationProblem(optf, collect(tunable_pas), prob_args)
     sol = Optimization.solve(optprob, opt.solve_alg, callback=callback_func, maxiters=opt.maxiters)
     opt_pas = update_ca(default_model_pas, ComponentVector(sol.u, tunable_axes))
@@ -133,7 +135,7 @@ function (opt::BatchOptimizer{C,S})(
     target::Vector;
     tunable_pas::ComponentVector,
     const_pas::ComponentVector,
-    config::Vector{<:NamedTuple}=fill(NamedTuple(), length(input)),
+    config::Vector=fill(NamedTuple(), length(input)),
     kwargs...
 ) where {C,S}
     loss_recorder = NamedTuple[]
@@ -154,6 +156,7 @@ function (opt::BatchOptimizer{C,S})(
 
     #* Constructing and solving optimization problems
     optf = Optimization.OptimizationFunction(opt.objective_func, opt.adtype)
+    @info "The size of tunable parameters is $(length(tunable_pas))"
     optprob = Optimization.OptimizationProblem(optf, collect(tunable_pas), prob_args)
     sol = Optimization.solve(optprob, opt.solve_alg, ncycle(train_batch, opt.maxiters), callback=callback_func)
     #* Returns the optimized model parameters
diff --git a/src/route.jl b/src/route.jl
index e743321..c30d04f 100644
--- a/src/route.jl
+++ b/src/route.jl
@@ -330,7 +330,7 @@ struct RapidRoute <: AbstractRapidRoute
         #* Extract all parameters names of funcs and dfuncs
         param_names = [:rapid_k, :rapid_x]
         #* Setup the name information of the hydrobucket
-        route_name = name === nothing ? Symbol(Symbol(reduce((x, y) -> Symbol(x, y), input_names)), :_vector_route) : name
+        route_name = name === nothing ? Symbol(Symbol(reduce((x, y) -> Symbol(x, y), input_names)), :_rapid_route) : name
         meta = HydroMeta(route_name, input_names, output_names, param_names, Symbol[], Symbol[])
         #* generate adjacency matrix from network
         adjacency = adjacency_matrix(network)'
diff --git a/src/solver.jl b/src/solver.jl
index 58731bb..324ac16 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -44,7 +44,7 @@ A custom ODEProblem solver
 """
 @kwdef struct DiscreteSolver <: AbstractHydroSolver
     alg = FunctionMap{true}()
-    sensealg = InterpolatingAdjoint()
+    sensealg = ZygoteAdjoint()
 end
 
 function (solver::DiscreteSolver)(
@@ -58,7 +58,7 @@ function (solver::DiscreteSolver)(
     #* build problem
     prob = DiscreteProblem(ode_func!, initstates, (timeidx[1], timeidx[end]), params)
     #* solve problem
-    sol = solve(prob, solver.alg, saveat=timeidx, sensealg=solver.sensealg)
+    sol = solve(prob, solver.alg, saveat=timeidx) # , sensealg=solver.sensealg
     if convert_to_array
         if SciMLBase.successful_retcode(sol)
             sol_arr = Array(sol)
@@ -90,8 +90,9 @@ struct ManualSolver{mutable} <: AbstractHydroSolver end
 function (solver::ManualSolver{true})(
     du_func::Function,
     pas::ComponentVector,
-    initstates::AbstractArray,
-    timeidx::AbstractVector
+    initstates::AbstractArray{<:Number, 1},
+    timeidx::AbstractVector;
+    convert_to_array::Bool=true
 )
     T1 = promote_type(eltype(pas), eltype(initstates))
     states_results = zeros(T1, size(initstates)..., length(timeidx))
@@ -104,18 +105,55 @@ function (solver::ManualSolver{true})(
     states_results
 end
 
+function (solver::ManualSolver{true})(
+    du_func::Function,
+    pas::ComponentVector,
+    initstates::AbstractArray{<:Number, 2},
+    timeidx::AbstractVector;
+    convert_to_array::Bool=true
+)
+    T1 = promote_type(eltype(pas), eltype(initstates))
+    states_results = zeros(T1, size(initstates)..., length(timeidx))
+    tmp_initstates = copy(initstates)
+    for (i, t) in enumerate(timeidx)
+        tmp_du = du_func(tmp_initstates, pas, t)
+        tmp_initstates = tmp_initstates .+ tmp_du
+        states_results[:, :, i] .= tmp_initstates
+    end
+    states_results
+end
+
 function (solver::ManualSolver{false})(
     du_func::Function,
     pas::ComponentVector,
-    initstates::AbstractArray,
-    timeidx::AbstractVector
+    initstates::AbstractArray{<:Number, 1},
+    timeidx::AbstractVector;
+    convert_to_array::Bool=true
 )
     states_results = []
     tmp_initstates = copy(initstates)
     for t in timeidx
         tmp_du = du_func(tmp_initstates, pas, t)
         tmp_initstates = tmp_initstates .+ tmp_du
-        states_results = vcat(states_results, tmp_initstates)
+        states_results = vcat(states_results, [tmp_initstates])
     end
     reduce((m1, m2) -> cat(m1, m2, dims=length(size(initstates))+1), states_results)
 end
+
+function (solver::ManualSolver{false})(
+    du_func::Function,
+    pas::ComponentVector,
+    initstates::AbstractArray{<:Number, 2},
+    timeidx::AbstractVector;
+    convert_to_array::Bool=true
+)
+    states_results = []
+    tmp_initstates = copy(initstates)
+    for t in timeidx
+        tmp_du = du_func(tmp_initstates, pas, t)
+        tmp_initstates = tmp_initstates .+ tmp_du
+        states_results = vcat(states_results, [tmp_initstates])
+    end
+    output = reduce((m1, m2) -> cat(m1, m2, dims=length(size(initstates))+1), states_results)
+    output
+end
diff --git a/src/uh.jl b/src/uh.jl
index d773b2b..3696830 100644
--- a/src/uh.jl
+++ b/src/uh.jl
@@ -102,7 +102,7 @@ struct UnitHydrograph{T<:Num,UF<:UHFunction,M<:HydroMeta,ST} <: AbstractHydrogra
         output::T,
         param::T;
         uhfunc::UF,
-        solvetype::Symbol=:DISCRETE,
+        solvetype::Symbol=:SPARSE,
     ) where {T<:Num,UF<:UHFunction}
         output_name = Symbolics.tosymbol(output, escape=false)
         @assert solvetype in [:DISCRETE, :SPARSE, :INTEGRAL] "solvetype must be one of [:DISCRETE, :SPARSE, :INTEGRAL]"
@@ -111,6 +111,17 @@ struct UnitHydrograph{T<:Num,UF<:UHFunction,M<:HydroMeta,ST} <: AbstractHydrogra
 
         return new{T,UF,typeof(meta),solvetype}([input], [output], [param], uhfunc, meta)
     end
+
+    function UnitHydrograph(
+        fluxes::Pair{Vector{T},Vector{T}},
+        param::T;
+        uhfunc::UF,
+        solvetype::Symbol=:SPARSE,
+    ) where {T<:Num,UF<:UHFunction}
+        input = fluxes[1][1]
+        output = fluxes[2][1]
+        return new{T,UF,typeof(meta),solvetype}(input, output, [param], uhfunc, meta)
+    end
 end
 
 """
diff --git a/src/utils/io.jl b/src/utils/io.jl
new file mode 100644
index 0000000..609ad51
--- /dev/null
+++ b/src/utils/io.jl
@@ -0,0 +1,38 @@
+struct NamedTupleIOAdapter{C,M} <: AbstractIOAdapter where {C<:AbstractComponent,M<:HydroMeta}
+    component::C
+    meta::M
+
+    function NamedTupleIOAdapter(component::C) where {C<:AbstractComponent}
+        new{C,typeof(component.meta)}(component, component.meta)
+    end
+end
+
+function (adapter::NamedTupleIOAdapter)(
+    input::NamedTuple,
+    pas::ComponentVector;
+    config::NamedTuple=NamedTuple(),
+    kwargs...
+)::NamedTuple
+    @assert((all(input_name in keys(input) for input_name in get_input_names(adapter))),
+        "Missing required inputs. Expected all of $(get_input_names(adapter)), but got $(keys(input)).")
+    input_matrix = Matrix(reduce(hcat, [input[k] for k in get_input_names(adapter)])')
+    output_matrix = adapter.component(input_matrix, pas; config=config, kwargs...)
+    return NamedTuple{Tuple(vcat(get_state_names(adapter), get_output_names(adapter)))}(eachslice(output_matrix, dims=1))
+end
+
+function (adapter::NamedTupleIOAdapter)(
+    input::Vector{<:NamedTuple},
+    pas::ComponentVector;
+    config::NamedTuple=NamedTuple(),
+    kwargs...
+)::Vector{<:NamedTuple}
+    for i in eachindex(input)
+        @assert(all(input_name in keys(input[i]) for input_name in get_input_names(adapter)),
+            "Missing required inputs. Expected all of $(get_input_names(adapter)), but got $(keys(input[i])) at $i input.")
+    end
+    input_mats = [reduce(hcat, collect(input[i][k] for k in get_input_names(adapter))) for i in eachindex(input)]
+    input_arr = reduce((m1, m2) -> cat(m1, m2, dims=3), input_mats)
+    output_arr = adapter.component(input_arr, pas; config=config, kwargs...)
+    output_names_tuple = Tuple(vcat(get_state_names(adapter), get_output_names(adapter)))
+    return [NamedTuple{output_names_tuple}(eachslice(output_arr_, dims=1)) for output_arr_ in eachslice(output_arr, dims=2)]
+end
\ No newline at end of file
diff --git a/src/utils/show.jl b/src/utils/show.jl
index 7e81b18..22d3814 100644
--- a/src/utils/show.jl
+++ b/src/utils/show.jl
@@ -3,17 +3,17 @@ function Base.show(io::IO, flux::AbstractHydroFlux)
 
     if compact
         print(io, "HydroFlux(")
-        print(io, "inputs: ", isempty(flux.meta.input) ? "nothing" : join(flux.meta.input, ", "))
-        print(io, ", outputs: ", isempty(flux.meta.output) ? "nothing" : join(flux.meta.output, ", "))
-        print(io, ", params: ", isempty(flux.meta.param) ? "nothing" : join(flux.meta.param, ", "))
+        print(io, "inputs: ", isempty(flux.meta.inputs) ? "nothing" : join(flux.meta.inputs, ", "))
+        print(io, ", outputs: ", isempty(flux.meta.outputs) ? "nothing" : join(flux.meta.outputs, ", "))
+        print(io, ", params: ", isempty(flux.meta.params) ? "nothing" : join(flux.meta.params, ", "))
         print(io, ")")
     else
         println(io, "HydroFlux:")
-        println(io, "  Inputs: ", isempty(flux.meta.input) ? "nothing" : join(flux.meta.input, ", "))
-        println(io, "  Outputs: ", isempty(flux.meta.output) ? "nothing" : join(flux.meta.output, ", "))
-        println(io, "  Parameters: ", isempty(flux.meta.param) ? "nothing" : join(flux.meta.param, ", "))
+        println(io, "  Inputs: ", isempty(flux.meta.inputs) ? "nothing" : join(flux.meta.inputs, ", "))
+        println(io, "  Outputs: ", isempty(flux.meta.outputs) ? "nothing" : join(flux.meta.outputs, ", "))
+        println(io, "  Parameters: ", isempty(flux.meta.params) ? "nothing" : join(flux.meta.params, ", "))
         println(io, "  Expressions:")
-        for (output, expr) in zip(flux.meta.output, flux.exprs)
+        for (output, expr) in zip(flux.meta.outputs, flux.exprs)
             println(io, "    $output = $expr")
         end
     end
@@ -25,15 +25,15 @@ function Base.show(io::IO, flux::AbstractStateFlux)
 
     if compact
         print(io, "StateFlux(")
-        print(io, "inputs: ", isempty(flux.meta.input) ? "nothing" : join(flux.meta.input, ", "))
-        print(io, ", params: ", isempty(flux.meta.param) ? "nothing" : join(flux.meta.param, ", "))
-        print(io, ", states: ", isempty(flux.meta.state) ? "nothing" : join(flux.meta.state, ", "))
+        print(io, "inputs: ", isempty(flux.meta.inputs) ? "nothing" : join(flux.meta.inputs, ", "))
+        print(io, ", params: ", isempty(flux.meta.params) ? "nothing" : join(flux.meta.params, ", "))
+        print(io, ", states: ", isempty(flux.meta.states) ? "nothing" : join(flux.meta.states, ", "))
         print(io, ")")
     else
         println(io, "StateFlux:")
-        println(io, "  Inputs: ", isempty(flux.meta.input) ? "nothing" : join(flux.meta.input, ", "))
-        println(io, "  Parameters: ", isempty(flux.meta.param) ? "nothing" : join(flux.meta.param, ", "))
-        println(io, "  States: ", isempty(flux.meta.state) ? "nothing" : join(flux.meta.state, ", "))
+        println(io, "  Inputs: ", isempty(flux.meta.inputs) ? "nothing" : join(flux.meta.inputs, ", "))
+        println(io, "  Parameters: ", isempty(flux.meta.params) ? "nothing" : join(flux.meta.params, ", "))
+        println(io, "  States: ", isempty(flux.meta.states) ? "nothing" : join(flux.meta.states, ", "))
         println(io, "  Expressions:")
         println(io, "    $(flux.state) = $(flux.expr)")
     end
@@ -142,29 +142,14 @@ function Base.show(io::IO, model::AbstractModel)
     buckets_in_model = filter(x -> x isa AbstractBucket, model.components)
     routes_in_model = filter(x -> x isa AbstractRoute, model.components)
     @assert(length(routes_in_model) <= 1, "Only one route is allowed in a model")
-    model_type = :LumpedModel
-    if length(routes_in_model) == 1
-        if routes_in_model[1] isa AbstractGridRoute
-            model_type = :GridModel
-        elseif routes_in_model[1] isa AbstractVectorRoute
-            model_type = :VectorModel
-        elseif routes_in_model[1] isa AbstractSumRoute
-            model_type = :MultiNodesModel
-        else
-            @error "Unknown route type: $(typeof(routes_in_model[1]))"
-        end
-    end
-
     compact = get(io, :compact, false)
     if compact
         print(io, "HydroModel(")
         print(io, "name: ", model.meta.name)
-        print(io, ", type: ", model_type)
         print(io, ", components: ", length(model.components))
         print(io, ")")
     else
         println(io, "HydroModel: ", model.meta.name)
-        println(io, "  Type: ", model_type)
         println(io, "  Components: ", join(map(c -> c.meta.name, model.components), ", "))
         println(io, "  Inputs: ", isempty(model.meta.inputs) ? "nothing" : join(model.meta.inputs, ", "))
         println(io, "  Outputs: ", isempty(model.meta.outputs) ? "nothing" : join(model.meta.outputs, ", "))
diff --git a/src/wrapper.jl b/src/wrapper.jl
index a1c36a1..f95f180 100644
--- a/src/wrapper.jl
+++ b/src/wrapper.jl
@@ -184,21 +184,18 @@ This is useful for focusing analysis on a particular outlet point in a hydrologi
 """
 struct ComputeComponentOutlet{C<:AbstractComponent,M<:HydroMeta} <: AbstractHydroWrapper
     component::C
-    outlet::Symbol
+    outlet::Integer
     meta::M
 
-    function ComputeComponentOutlet(component::C, outlet::S) where {C<:AbstractComponent,S<:Symbol}
-        @assert outlet in component.hrunames "Outlet $outlet is not a valid output name for component $component"
-        new{C,typeof(meta)}(component, outlet, meta)
+    function ComputeComponentOutlet(component::C, outlet::I) where {C<:AbstractComponent,I<:Integer}
+        new{C,typeof(component.meta)}(component, outlet, component.meta)
     end
 end
 
 function (wrapper::ComputeComponentOutlet)(input::AbstractVector{<:NamedTuple}, pas::ComponentVector; kwargs...)
     convert_to_ntp = get(kwargs, :convert_to_ntp, true)
-    @assert wrapper.outlet in wrapper.component.varnames "Outlet $wrapper.outlet is not a valid output name for component $wrapper.component"
-    outlet_idx = findfirst(nm -> nm == wrapper.outlet, wrapper.component.varnames)
     output = wrapper.component(input, pas; kwargs...)
-    return convert_to_ntp ? output[outlet_idx] : output[:, outlet_idx, :]
+    return convert_to_ntp ? output[wrapper.outlet] : output[:, wrapper.outlet, :]
 end
 
 """
@@ -285,11 +282,11 @@ end
 
 function (wrapper::NeuralWrapper)(input::Array{T,3}, pas::ComponentVector; kwargs...) where {T}
     nn_ps = pas[:nn][wrapper.meta.name]
-    output = wrapper.func.(eachslice(input, dims=2), Ref(nn_ps))
+    output = wrapper.func.(Matrix.(eachslice(input, dims=2)), Ref(nn_ps))
     convert_to_ntp = get(kwargs, :convert_to_ntp, false)
     if convert_to_ntp
-        return [NamedTuple{Tuple(get_output_names(wrapper))}(eachslice(output[:, i, :], dims=1)) for i in axes(output, 2)]
+        return [NamedTuple{Tuple(get_output_names(wrapper))}(eachslice(output_, dims=1)) for output_ in eachslice(output, dims=2)]
     else
-        return output
+        return permutedims(reduce((m1, m2) -> cat(m1, m2, dims=3), output), (1,3,2))
     end
 end
\ No newline at end of file