diff --git a/necsim/impls/no-std/src/array2d.rs b/necsim/impls/no-std/src/array2d.rs
index dd4552ebe..e32a2a14e 100644
--- a/necsim/impls/no-std/src/array2d.rs
+++ b/necsim/impls/no-std/src/array2d.rs
@@ -4,12 +4,17 @@
 /// Based on a subset of Harrison McCullough's MIT-licensed [`array2d`] crate.
 ///
 /// [`array2d`]: https://github.com/HarrisonMc555/array2d
-use alloc::{boxed::Box, vec::Vec};
+use alloc::{boxed::Box, rc::Rc, sync::Arc, vec::Vec};
 
-use core::ops::{Index, IndexMut};
+use core::{
+    marker::PhantomData,
+    ops::{Deref, DerefMut, Index, IndexMut},
+};
+
+pub trait ArrayBackend<T>: From<Vec<T>> + Deref<Target = [T]> {}
+impl<T, B: From<Vec<T>> + Deref<Target = [T]>> ArrayBackend<T> for B {}
 
 /// A fixed sized two-dimensional array.
-#[derive(Clone, Eq, PartialEq)]
 #[cfg_attr(feature = "cuda", derive(rust_cuda::lend::LendRustToCuda))]
 #[cfg_attr(
     feature = "cuda",
@@ -18,13 +23,36 @@ use core::ops::{Index, IndexMut};
         bound = "T: rust_cuda::safety::PortableBitSemantics + const_type_layout::TypeGraphLayout"
     )
 )]
-pub struct Array2D<T> {
+pub struct Array2D<T, B: ArrayBackend<T> = Box<[T]>> {
     #[cfg_attr(feature = "cuda", cuda(embed))]
-    array: Box<[T]>,
+    array: B,
     num_rows: usize,
+    marker: PhantomData<T>,
+}
+
+impl<T, B: ArrayBackend<T> + Clone> Clone for Array2D<T, B> {
+    fn clone(&self) -> Self {
+        Self {
+            array: self.array.clone(),
+            num_rows: self.num_rows,
+            marker: PhantomData::<T>,
+        }
+    }
+}
+
+impl<T, B: ArrayBackend<T> + PartialEq> PartialEq for Array2D<T, B> {
+    fn eq(&self, other: &Self) -> bool {
+        (self.array == other.array) && (self.num_rows == other.num_rows)
+    }
 }
+impl<T, B: ArrayBackend<T> + Eq> Eq for Array2D<T, B> {}
+
+pub type ArcArray2D<T> = Array2D<T, Arc<[T]>>;
+pub type BoxArray2D<T> = Array2D<T, Box<[T]>>;
+pub type RcArray2D<T> = Array2D<T, Rc<[T]>>;
+pub type VecArray2D<T> = Array2D<T, Vec<T>>;
 
-impl<T> core::fmt::Debug for Array2D<T> {
+impl<T, B: ArrayBackend<T>> core::fmt::Debug for Array2D<T, B> {
     fn fmt(&self, fmt: &mut core::fmt::Formatter) -> core::fmt::Result {
         write!(
             fmt,
@@ -49,7 +77,7 @@ pub enum Error {
     NotEnoughElements,
 }
 
-impl<T> Array2D<T> {
+impl<T, B: ArrayBackend<T>> Array2D<T, B> {
     /// Creates a new [`Array2D`] from a slice of rows, each of which is a
     /// [`Vec`] of elements.
     ///
@@ -86,8 +114,9 @@ impl<T> Array2D<T> {
                 .iter()
                 .flat_map(Vec::clone)
                 .collect::<Vec<_>>()
-                .into_boxed_slice(),
+                .into(),
             num_rows: elements.len(),
+            marker: PhantomData::<T>,
         })
     }
 
@@ -131,8 +160,9 @@ impl<T> Array2D<T> {
             return Err(Error::DimensionMismatch);
         }
         Ok(Array2D {
-            array: elements.to_vec().into_boxed_slice(),
+            array: elements.to_vec().into(),
             num_rows,
+            marker: PhantomData::<T>,
         })
     }
 
@@ -155,8 +185,9 @@ impl<T> Array2D<T> {
         let total_len = num_rows * num_columns;
         let array = alloc::vec![element; total_len];
         Array2D {
-            array: array.into_boxed_slice(),
+            array: array.into(),
             num_rows,
+            marker: PhantomData::<T>,
         }
     }
 
@@ -201,8 +232,9 @@ impl<T> Array2D<T> {
             return Err(Error::NotEnoughElements);
         }
         Ok(Array2D {
-            array: array.into_boxed_slice(),
+            array: array.into(),
             num_rows,
+            marker: PhantomData::<T>,
         })
     }
 
@@ -279,7 +311,10 @@ impl<T> Array2D<T> {
     ///
     /// [`Some`]: https://doc.rust-lang.org/std/option/enum.Option.html#variant.Some
     /// [`None`]: https://doc.rust-lang.org/std/option/enum.Option.html#variant.None
-    pub fn get_mut(&mut self, row: usize, column: usize) -> Option<&mut T> {
+    pub fn get_mut(&mut self, row: usize, column: usize) -> Option<&mut T>
+    where
+        B: DerefMut,
+    {
         self.get_index(row, column)
             .map(move |index| &mut self.array[index])
     }
@@ -412,10 +447,62 @@ impl<T> Array2D<T> {
     /// [`Vec`]: https://doc.rust-lang.org/std/vec/struct.Vec.html
     /// [row major order]: https://en.wikipedia.org/wiki/Row-_and_column-major_order
     #[must_use]
-    pub fn into_row_major(self) -> Vec<T> {
+    pub fn into_row_major(self) -> Vec<T>
+    where
+        B: Into<Vec<T>>,
+    {
         self.array.into()
     }
 
+    /// Converts the [`Array2D`] into its inner [`ArrayBackend`] `B` of
+    /// elements in [row major order].
+    ///
+    /// # Examples
+    ///
+    /// ```
+    /// # use necsim_impls_no_std::array2d::{VecArray2D, Error};
+    /// # fn main() -> Result<(), Error> {
+    /// let rows = vec![vec![1, 2, 3], vec![4, 5, 6]];
+    /// let array = VecArray2D::from_rows(&rows)?;
+    /// assert_eq!(array.into_row_major_inner(), vec![1, 2, 3, 4, 5, 6]);
+    /// # Ok(())
+    /// # }
+    /// ```
+    ///
+    /// [`Array2D`]: struct.Array2D.html
+    /// [row major order]: https://en.wikipedia.org/wiki/Row-_and_column-major_order
+    #[must_use]
+    pub fn into_row_major_inner(self) -> B {
+        self.array
+    }
+
+    /// Converts the [`Array2D`] into from its current into a new
+    /// [`ArrayBackend`].
+    ///
+    /// # Examples
+    ///
+    /// ```
+    /// # use necsim_impls_no_std::array2d::{VecArray2D, Error};
+    /// # fn main() -> Result<(), Error> {
+    /// let rows = vec![vec![1, 2, 3], vec![4, 5, 6]];
+    /// let array = VecArray2D::from_rows(&rows)?;
+    /// let array: BoxArray2D = array.into_backend();
+    /// assert_eq!(array.into_row_major(), vec![1, 2, 3, 4, 5, 6]);
+    /// # Ok(())
+    /// # }
+    /// ```
+    ///
+    /// [`Array2D`]: struct.Array2D.html
+    /// [row major order]: https://en.wikipedia.org/wiki/Row-_and_column-major_order
+    #[must_use]
+    pub fn switch_backend<B2: ArrayBackend<T> + From<B>>(self) -> Array2D<T, B2> {
+        Array2D {
+            array: self.array.into(),
+            num_rows: self.num_rows,
+            marker: PhantomData::<T>,
+        }
+    }
+
     /// Returns an [`Iterator`] over references to all elements in
     /// [row major order].
     ///
@@ -449,7 +536,7 @@ impl<T> Array2D<T> {
     }
 }
 
-impl<T> Index<(usize, usize)> for Array2D<T> {
+impl<T, B: ArrayBackend<T>> Index<(usize, usize)> for Array2D<T, B> {
     type Output = T;
 
     /// Returns the element at the given indices, given as `(row, column)`.
@@ -477,7 +564,7 @@ impl<T> Index<(usize, usize)> for Array2D<T> {
     }
 }
 
-impl<T> IndexMut<(usize, usize)> for Array2D<T> {
+impl<T, B: ArrayBackend<T> + DerefMut> IndexMut<(usize, usize)> for Array2D<T, B> {
     /// Returns a mutable version of the element at the given indices, given as
     /// `(row, column)`.
     ///
diff --git a/necsim/impls/no-std/src/cogs/dispersal_sampler/in_memory/alias/mod.rs b/necsim/impls/no-std/src/cogs/dispersal_sampler/in_memory/alias/mod.rs
index 5b2c5a823..d2db471e2 100644
--- a/necsim/impls/no-std/src/cogs/dispersal_sampler/in_memory/alias/mod.rs
+++ b/necsim/impls/no-std/src/cogs/dispersal_sampler/in_memory/alias/mod.rs
@@ -9,7 +9,8 @@ use necsim_core::{
 use necsim_core_bond::NonNegativeF64;
 
 use crate::{
-    alias::AliasMethodSampler, array2d::Array2D,
+    alias::AliasMethodSampler,
+    array2d::{ArcArray2D, Array2D},
     cogs::dispersal_sampler::in_memory::InMemoryDispersalSampler,
 };
 
@@ -18,8 +19,7 @@ mod dispersal;
 #[allow(clippy::module_name_repetitions)]
 #[derive(Debug)]
 pub struct InMemoryAliasDispersalSampler<M: MathsCore, H: Habitat<M>, G: RngCore<M>> {
-    // TODO: use Arc
-    alias_dispersal: Array2D<Option<AliasMethodSampler<usize>>>,
+    alias_dispersal: ArcArray2D<Option<AliasMethodSampler<usize>>>,
     marker: PhantomData<(M, H, G)>,
 }
 
diff --git a/necsim/impls/no-std/src/cogs/dispersal_sampler/in_memory/packed_alias/mod.rs b/necsim/impls/no-std/src/cogs/dispersal_sampler/in_memory/packed_alias/mod.rs
index 08bc4a414..99a96cf9e 100644
--- a/necsim/impls/no-std/src/cogs/dispersal_sampler/in_memory/packed_alias/mod.rs
+++ b/necsim/impls/no-std/src/cogs/dispersal_sampler/in_memory/packed_alias/mod.rs
@@ -2,14 +2,15 @@ use alloc::{sync::Arc, vec::Vec};
 use core::{marker::PhantomData, ops::Range};
 use necsim_core_bond::NonNegativeF64;
 
-use r#final::Final;
-
 use necsim_core::{
     cogs::{Habitat, MathsCore, RngCore},
     landscape::Location,
 };
 
-use crate::{alias::packed::AliasMethodSamplerAtom, array2d::Array2D};
+use crate::{
+    alias::packed::AliasMethodSamplerAtom,
+    array2d::{ArcArray2D, Array2D},
+};
 
 mod dispersal;
 
@@ -43,9 +44,8 @@ impl From<AliasSamplerRange> for Range<usize> {
 #[cfg_attr(feature = "cuda", derive(rust_cuda::lend::LendRustToCuda))]
 #[cfg_attr(feature = "cuda", cuda(free = "M", free = "H", free = "G"))]
 pub struct InMemoryPackedAliasDispersalSampler<M: MathsCore, H: Habitat<M>, G: RngCore<M>> {
-    // TODO: use Arc
     #[cfg_attr(feature = "cuda", cuda(embed))]
-    alias_dispersal_ranges: Final<Array2D<AliasSamplerRange>>,
+    alias_dispersal_ranges: ArcArray2D<AliasSamplerRange>,
     #[cfg_attr(feature = "cuda", cuda(embed))]
     alias_dispersal_buffer: Arc<[AliasMethodSamplerAtom<usize>]>,
     marker: PhantomData<(M, H, G)>,
@@ -107,7 +107,7 @@ impl<M: MathsCore, H: Habitat<M>, G: RngCore<M>> InMemoryDispersalSampler<M, H,
         .unwrap(); // infallible by PRE;
 
         Self {
-            alias_dispersal_ranges: Final::new(alias_dispersal_ranges),
+            alias_dispersal_ranges,
             alias_dispersal_buffer: Arc::from(alias_dispersal_buffer.into_boxed_slice()),
             marker: PhantomData::<(M, H, G)>,
         }
@@ -137,7 +137,7 @@ impl<M: MathsCore, H: Habitat<M>, G: RngCore<M>> Clone
 {
     fn clone(&self) -> Self {
         Self {
-            alias_dispersal_ranges: Final::new(self.alias_dispersal_ranges.clone()),
+            alias_dispersal_ranges: self.alias_dispersal_ranges.clone(),
             alias_dispersal_buffer: self.alias_dispersal_buffer.clone(),
             marker: PhantomData::<(M, H, G)>,
         }
diff --git a/necsim/impls/no-std/src/cogs/dispersal_sampler/in_memory/packed_separable_alias/mod.rs b/necsim/impls/no-std/src/cogs/dispersal_sampler/in_memory/packed_separable_alias/mod.rs
index a317da739..f65929e1a 100644
--- a/necsim/impls/no-std/src/cogs/dispersal_sampler/in_memory/packed_separable_alias/mod.rs
+++ b/necsim/impls/no-std/src/cogs/dispersal_sampler/in_memory/packed_separable_alias/mod.rs
@@ -2,14 +2,15 @@ use alloc::{sync::Arc, vec::Vec};
 use core::marker::PhantomData;
 use necsim_core_bond::{ClosedUnitF64, NonNegativeF64};
 
-use r#final::Final;
-
 use necsim_core::{
     cogs::{Habitat, MathsCore, RngCore},
     landscape::Location,
 };
 
-use crate::{alias::packed::AliasMethodSamplerAtom, array2d::Array2D};
+use crate::{
+    alias::packed::AliasMethodSamplerAtom,
+    array2d::{ArcArray2D, Array2D, VecArray2D},
+};
 
 mod dispersal;
 
@@ -42,12 +43,10 @@ pub struct SeparableAliasSelfDispersal {
 #[cfg_attr(feature = "cuda", cuda(free = "M", free = "H", free = "G"))]
 pub struct InMemoryPackedSeparableAliasDispersalSampler<M: MathsCore, H: Habitat<M>, G: RngCore<M>>
 {
-    // TODO: use Arc
     #[cfg_attr(feature = "cuda", cuda(embed))]
-    alias_dispersal_ranges: Final<Array2D<AliasSamplerRange>>,
-    // TODO: use Arc
+    alias_dispersal_ranges: ArcArray2D<AliasSamplerRange>,
     #[cfg_attr(feature = "cuda", cuda(embed))]
-    self_dispersal: Final<Array2D<SeparableAliasSelfDispersal>>,
+    self_dispersal: ArcArray2D<SeparableAliasSelfDispersal>,
     #[cfg_attr(feature = "cuda", cuda(embed))]
     alias_dispersal_buffer: Arc<[AliasMethodSamplerAtom<usize>]>,
     marker: PhantomData<(M, H, G)>,
@@ -67,7 +66,7 @@ impl<M: MathsCore, H: Habitat<M>, G: RngCore<M>> InMemoryDispersalSampler<M, H,
 
         let mut alias_dispersal_buffer = Vec::new();
 
-        let mut self_dispersal = Array2D::filled_with(
+        let mut self_dispersal = VecArray2D::filled_with(
             SeparableAliasSelfDispersal {
                 self_dispersal: ClosedUnitF64::zero(),
                 non_self_dispersal_event: usize::MAX,
@@ -180,8 +179,8 @@ impl<M: MathsCore, H: Habitat<M>, G: RngCore<M>> InMemoryDispersalSampler<M, H,
         .unwrap(); // infallible by PRE;
 
         Self {
-            alias_dispersal_ranges: Final::new(alias_dispersal_ranges),
-            self_dispersal: Final::new(self_dispersal),
+            alias_dispersal_ranges,
+            self_dispersal: self_dispersal.switch_backend(),
             alias_dispersal_buffer: Arc::from(alias_dispersal_buffer.into_boxed_slice()),
             marker: PhantomData::<(M, H, G)>,
         }
@@ -202,8 +201,8 @@ impl<M: MathsCore, H: Habitat<M>, G: RngCore<M>> Clone
 {
     fn clone(&self) -> Self {
         Self {
-            alias_dispersal_ranges: Final::new(self.alias_dispersal_ranges.clone()),
-            self_dispersal: Final::new(self.self_dispersal.clone()),
+            alias_dispersal_ranges: self.alias_dispersal_ranges.clone(),
+            self_dispersal: self.self_dispersal.clone(),
             alias_dispersal_buffer: self.alias_dispersal_buffer.clone(),
             marker: PhantomData::<(M, H, G)>,
         }
diff --git a/necsim/impls/no-std/src/cogs/dispersal_sampler/in_memory/separable_alias/mod.rs b/necsim/impls/no-std/src/cogs/dispersal_sampler/in_memory/separable_alias/mod.rs
index 7a90a87ca..27338a699 100644
--- a/necsim/impls/no-std/src/cogs/dispersal_sampler/in_memory/separable_alias/mod.rs
+++ b/necsim/impls/no-std/src/cogs/dispersal_sampler/in_memory/separable_alias/mod.rs
@@ -9,7 +9,8 @@ use necsim_core::{
 use necsim_core_bond::{ClosedUnitF64, NonNegativeF64};
 
 use crate::{
-    alias::AliasMethodSampler, array2d::Array2D,
+    alias::AliasMethodSampler,
+    array2d::{ArcArray2D, Array2D, VecArray2D},
     cogs::dispersal_sampler::in_memory::InMemoryDispersalSampler,
 };
 
@@ -18,10 +19,8 @@ mod dispersal;
 #[allow(clippy::module_name_repetitions)]
 #[derive(Debug)]
 pub struct InMemorySeparableAliasDispersalSampler<M: MathsCore, H: Habitat<M>, G: RngCore<M>> {
-    // TODO: use Arc
-    alias_dispersal: Array2D<Option<AliasMethodSampler<usize>>>,
-    // TODO: use Arc
-    self_dispersal: Array2D<ClosedUnitF64>,
+    alias_dispersal: ArcArray2D<Option<AliasMethodSampler<usize>>>,
+    self_dispersal: ArcArray2D<ClosedUnitF64>,
     _marker: PhantomData<(M, H, G)>,
 }
 
@@ -37,7 +36,7 @@ impl<M: MathsCore, H: Habitat<M>, G: RngCore<M>> InMemoryDispersalSampler<M, H,
         let mut event_weights: Vec<(usize, NonNegativeF64)> =
             Vec::with_capacity(dispersal.row_len());
 
-        let mut self_dispersal = Array2D::filled_with(
+        let mut self_dispersal = VecArray2D::filled_with(
             ClosedUnitF64::zero(),
             usize::from(habitat_extent.height()),
             usize::from(habitat_extent.width()),
@@ -108,7 +107,7 @@ impl<M: MathsCore, H: Habitat<M>, G: RngCore<M>> InMemoryDispersalSampler<M, H,
 
         Self {
             alias_dispersal,
-            self_dispersal,
+            self_dispersal: self_dispersal.switch_backend(),
             _marker: PhantomData::<(M, H, G)>,
         }
     }