diff --git a/docs/source/api/exports.rst b/docs/source/api/exports.rst
index a064f1fd6..92ee969c4 100644
--- a/docs/source/api/exports.rst
+++ b/docs/source/api/exports.rst
@@ -47,6 +47,14 @@ Exports
     :members:
     :show-inheritance:
 
+.. autoclass:: SurfaceFluxCylindrical
+    :members:
+    :show-inheritance:   
+
+.. autoclass:: SurfaceFluxSpherical
+    :members:
+    :show-inheritance:    
+
 .. autoclass:: HydrogenFlux
     :members:
     :show-inheritance:
diff --git a/festim/__init__.py b/festim/__init__.py
index 91118cb9e..bd29c6078 100644
--- a/festim/__init__.py
+++ b/festim/__init__.py
@@ -67,7 +67,11 @@
     VolumeQuantity,
     SurfaceQuantity,
 )
-from .exports.derived_quantities.surface_flux import SurfaceFlux
+from .exports.derived_quantities.surface_flux import (
+    SurfaceFlux,
+    SurfaceFluxCylindrical,
+    SurfaceFluxSpherical,
+)
 from .exports.derived_quantities.hydrogen_flux import HydrogenFlux
 from .exports.derived_quantities.thermal_flux import ThermalFlux
 from .exports.derived_quantities.average_volume import AverageVolume
diff --git a/festim/exports/derived_quantities/surface_flux.py b/festim/exports/derived_quantities/surface_flux.py
index a5f0786c3..f4ea91daa 100644
--- a/festim/exports/derived_quantities/surface_flux.py
+++ b/festim/exports/derived_quantities/surface_flux.py
@@ -1,5 +1,6 @@
 from festim import SurfaceQuantity, k_B
 import fenics as f
+import numpy as np
 
 
 class SurfaceFlux(SurfaceQuantity):
@@ -22,14 +23,17 @@ def __init__(self, field, surface) -> None:
         super().__init__(field=field, surface=surface)
         self.title = "Flux surface {}: {}".format(self.surface, self.field)
 
-    def compute(self, soret=False):
+    @property
+    def prop(self):
         field_to_prop = {
             "0": self.D,
             "solute": self.D,
             0: self.D,
             "T": self.thermal_cond,
         }
-        self.prop = field_to_prop[self.field]
+        return field_to_prop[self.field]
+
+    def compute(self, soret=False):
         flux = f.assemble(
             self.prop * f.dot(f.grad(self.function), self.n) * self.ds(self.surface)
         )
@@ -43,3 +47,141 @@ def compute(self, soret=False):
                 * self.ds(self.surface)
             )
         return flux
+
+
+class SurfaceFluxCylindrical(SurfaceFlux):
+    """
+    Object to compute the flux J of a field u through a surface
+    J = integral(-prop * grad(u) . n ds)
+    where prop is the property of the field (D, thermal conductivity, etc)
+    u is the field
+    n is the normal vector of the surface
+    ds is the surface measure in cylindrical coordinates.
+    ds = r dr dtheta or ds = r dz dtheta
+
+    Note: for particle fluxes J is given in H/s, for heat fluxes J is given in W
+
+    Args:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+        azimuth_range (tuple, optional): Range of the azimuthal angle
+            (theta) needs to be between 0 and 2 pi. Defaults to (0, 2 * np.pi).
+    """
+
+    def __init__(self, field, surface, azimuth_range=(0, 2 * np.pi)) -> None:
+        super().__init__(field, surface)
+        self.r = None
+        self.azimuth_range = azimuth_range
+
+    @property
+    def azimuth_range(self):
+        return self._azimuth_range
+
+    @azimuth_range.setter
+    def azimuth_range(self, value):
+        if value[0] < 0 or value[1] > 2 * np.pi:
+            raise ValueError("Azimuthal range must be between 0 and pi")
+        self._azimuth_range = value
+
+    def compute(self, soret=False):
+        if soret:
+            raise NotImplementedError(
+                "Soret effect not implemented for cylindrical coordinates"
+            )
+
+        if self.r is None:
+            mesh = (
+                self.function.function_space().mesh()
+            )  # get the mesh from the function
+            rthetaz = f.SpatialCoordinate(mesh)  # get the coordinates from the mesh
+            self.r = rthetaz[0]  # only care about r here
+
+        # dS_z = r dr dtheta , assuming axisymmetry dS_z = theta r dr
+        # dS_r = r dz dtheta , assuming axisymmetry dS_r = theta r dz
+        # in both cases the expression with self.ds is the same
+        # we assume full cylinder theta = 2 pi
+        flux = f.assemble(
+            self.prop
+            * self.r
+            * f.dot(f.grad(self.function), self.n)
+            * self.ds(self.surface)
+        )
+        flux *= self.azimuth_range[1] - self.azimuth_range[0]
+        return flux
+
+
+class SurfaceFluxSpherical(SurfaceFlux):
+    """
+    Object to compute the flux J of a field u through a surface
+    J = integral(-prop * grad(u) . n ds)
+    where prop is the property of the field (D, thermal conductivity, etc)
+    u is the field
+    n is the normal vector of the surface
+    ds is the surface measure in spherical coordinates.
+    ds = r^2 sin(theta) dtheta dphi
+
+    Note: for particle fluxes J is given in H/s, for heat fluxes J is given in W
+
+    Args:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+        azimuth_range (tuple, optional): Range of the azimuthal angle
+            (phi) needs to be between 0 and pi. Defaults to (0, np.pi).
+        polar_range (tuple, optional): Range of the polar angle
+            (theta) needs to be between - pi and pi. Defaults to (-np.pi, np.pi).
+    """
+
+    def __init__(
+        self, field, surface, azimuth_range=(0, np.pi), polar_range=(-np.pi, np.pi)
+    ) -> None:
+        super().__init__(field, surface)
+        self.r = None
+        self.polar_range = polar_range
+        self.azimuth_range = azimuth_range
+
+    @property
+    def polar_range(self):
+        return self._polar_range
+
+    @polar_range.setter
+    def polar_range(self, value):
+        if value[0] < -np.pi or value[1] > np.pi:
+            raise ValueError("Polar range must be between - pi and pi")
+        self._polar_range = value
+
+    @property
+    def azimuth_range(self):
+        return self._azimuth_range
+
+    @azimuth_range.setter
+    def azimuth_range(self, value):
+        if value[0] < 0 or value[1] > np.pi:
+            raise ValueError("Azimuthal range must be between 0 and pi")
+        self._azimuth_range = value
+
+    def compute(self, soret=False):
+        if soret:
+            raise NotImplementedError(
+                "Soret effect not implemented for spherical coordinates"
+            )
+
+        if self.r is None:
+            mesh = (
+                self.function.function_space().mesh()
+            )  # get the mesh from the function
+            rthetaphi = f.SpatialCoordinate(mesh)  # get the coordinates from the mesh
+            self.r = rthetaphi[0]  # only care about r here
+
+        # dS_r = r^2 sin(theta) dtheta dphi
+        # integral(f dS_r) = integral(f r^2 sin(theta) dtheta dphi)
+        #                  = (phi2 - phi1) * (-cos(theta2) + cos(theta1)) * f r^2
+        flux = f.assemble(
+            self.prop
+            * self.r**2
+            * f.dot(f.grad(self.function), self.n)
+            * self.ds(self.surface)
+        )
+        flux *= (self.polar_range[1] - self.polar_range[0]) * (
+            -np.cos(self.azimuth_range[1]) + np.cos(self.azimuth_range[0])
+        )
+        return flux
diff --git a/festim/generic_simulation.py b/festim/generic_simulation.py
index fe149271f..e61db7b54 100644
--- a/festim/generic_simulation.py
+++ b/festim/generic_simulation.py
@@ -297,25 +297,32 @@ def initialise(self):
 
         self.h_transport_problem.initialise(self.mesh, self.materials, self.dt)
 
-        # raise warning if SurfaceFlux is used with non-cartesian meshes
+        # raise warning if the derived quantities don't match the type of mesh
+        # eg. SurfaceFlux is used with cylindrical mesh
+        all_types_quantities = [
+            festim.MaximumSurface,
+            festim.MinimumSurface,
+            festim.MaximumVolume,
+            festim.MinimumVolume,
+            festim.PointValue,
+        ]  # these quantities can be used with any mesh
+        allowed_quantities = {
+            "cartesian": [
+                festim.SurfaceFlux,
+                festim.AverageSurface,
+                festim.AverageVolume,
+            ]
+            + all_types_quantities,
+            "cylindrical": [festim.SurfaceFluxCylindrical] + all_types_quantities,
+            "spherical": [festim.SurfaceFluxSpherical] + all_types_quantities,
+        }
+
         for export in self.exports:
             if isinstance(export, festim.DerivedQuantities):
-                allowed_quantities = (
-                    festim.MaximumSurface,
-                    festim.MinimumSurface,
-                    festim.MaximumVolume,
-                    festim.MinimumVolume,
-                    festim.PointValue,
-                )
-                if any(
-                    [
-                        not isinstance(q, allowed_quantities)
-                        for q in export.derived_quantities
-                    ]
-                ):
-                    if self.mesh.type != "cartesian":
+                for q in export.derived_quantities:
+                    if not isinstance(q, tuple(allowed_quantities[self.mesh.type])):
                         warnings.warn(
-                            "Some derived quantities may not work as intended for non-cartesian meshes"
+                            f"{type(q)} may not work as intended for {self.mesh.type} meshes"
                         )
         self.exports.initialise_derived_quantities(
             self.mesh.dx, self.mesh.ds, self.materials
diff --git a/test/simulation/test_initialise.py b/test/simulation/test_initialise.py
index 6f6ad95bc..da78696f0 100644
--- a/test/simulation/test_initialise.py
+++ b/test/simulation/test_initialise.py
@@ -137,5 +137,5 @@ def test_cartesian_and_surface_flux_warning(quantity, sys):
     my_model.exports = [derived_quantities]
 
     # test
-    with pytest.warns(UserWarning, match="Some derived quantities .* non-cartesian"):
+    with pytest.warns(UserWarning, match=f"may not work as intended for {sys} meshes"):
         my_model.initialise()
diff --git a/test/unit/test_exports/test_derived_quantities/test_surface_flux.py b/test/unit/test_exports/test_derived_quantities/test_surface_flux.py
index 3397049b5..668d805e6 100644
--- a/test/unit/test_exports/test_derived_quantities/test_surface_flux.py
+++ b/test/unit/test_exports/test_derived_quantities/test_surface_flux.py
@@ -1,5 +1,8 @@
-from festim import SurfaceFlux, k_B
+from festim import SurfaceFlux, k_B, SurfaceFluxCylindrical, SurfaceFluxSpherical
 import fenics as f
+import math
+import numpy as np
+import pytest
 
 
 def test_title_H():
@@ -80,3 +83,181 @@ def test_h_flux_with_soret(self):
         )
         flux = self.my_h_flux.compute(soret=True)
         assert flux == expected_flux
+
+
+@pytest.mark.parametrize("radius", [2, 3])
+@pytest.mark.parametrize("r0", [0, 2])
+@pytest.mark.parametrize("height", [2, 3])
+@pytest.mark.parametrize("c_top", [2, 3])
+@pytest.mark.parametrize("c_bottom", [2, 3])
+def test_compute_cylindrical(r0, radius, height, c_top, c_bottom):
+    """
+    Test that SurfaceFluxCylindrical computes the flux correctly on a hollow cylinder
+
+    Args:
+        r0 (float): internal radius
+        radius (float): cylinder radius
+        height (float): cylinder height
+        c_top (float): concentration top
+        c_bottom (float): concentration bottom
+    """
+    # creating a mesh with FEniCS
+    r1 = r0 + radius
+    z0 = 0
+    z1 = z0 + height
+    mesh_fenics = f.RectangleMesh(f.Point(r0, z0), f.Point(r1, z1), 10, 10)
+
+    # marking physical groups (volumes and surfaces)
+    volume_markers = f.MeshFunction("size_t", mesh_fenics, mesh_fenics.topology().dim())
+    volume_markers.set_all(1)
+
+    bot_surface = f.CompiledSubDomain(
+        f"on_boundary && near(x[1], {z0}, tol)", tol=1e-14
+    )
+    top_surface = f.CompiledSubDomain(
+        f"on_boundary && near(x[1], {z1}, tol)", tol=1e-14
+    )
+
+    surface_markers = f.MeshFunction(
+        "size_t", mesh_fenics, mesh_fenics.topology().dim() - 1
+    )
+    surface_markers.set_all(0)
+    ds = f.Measure("ds", domain=mesh_fenics, subdomain_data=surface_markers)
+    # Surface ids
+    bottom_id = 1
+    top_id = 2
+    bot_surface.mark(surface_markers, bottom_id)
+    top_surface.mark(surface_markers, top_id)
+
+    my_flux = SurfaceFluxCylindrical("solute", top_id)
+    my_flux.D = f.Constant(2)
+    V = f.FunctionSpace(mesh_fenics, "P", 1)
+    expr = f.Expression(
+        f"{c_bottom} + {(c_top - c_bottom)/(z1-z0)} * x[1]",
+        degree=1,
+    )
+    my_flux.function = f.interpolate(expr, V)
+
+    my_flux.n = f.FacetNormal(mesh_fenics)
+    my_flux.ds = ds
+
+    computed_value = my_flux.compute()
+    expected_value = (
+        -2
+        * math.pi
+        * float(my_flux.D)
+        * (c_bottom - c_top)
+        / height
+        * (0.5 * r1**2 - 0.5 * r0**2)
+    )
+    assert np.isclose(computed_value, expected_value)
+
+
+@pytest.mark.parametrize("r0", [1, 2])
+@pytest.mark.parametrize("radius", [1, 2])
+@pytest.mark.parametrize("c_left", [3, 4])
+@pytest.mark.parametrize("c_right", [1, 2])
+def test_compute_spherical(r0, radius, c_left, c_right):
+    """
+    Test that SurfaceFluxSpherical computes the flux correctly on a hollow sphere
+
+    Args:
+        r0 (float): internal radius
+        radius (float): sphere radius
+        c_left (float): concentration left
+        c_right (float): concentration right
+    """
+    # creating a mesh with FEniCS
+    r1 = r0 + radius
+    mesh_fenics = f.IntervalMesh(10, r0, r1)
+
+    # marking physical groups (volumes and surfaces)
+    volume_markers = f.MeshFunction("size_t", mesh_fenics, mesh_fenics.topology().dim())
+    volume_markers.set_all(1)
+
+    left_surface = f.CompiledSubDomain(
+        f"on_boundary && near(x[0], {r0}, tol)", tol=1e-14
+    )
+    right_surface = f.CompiledSubDomain(
+        f"on_boundary && near(x[0], {r1}, tol)", tol=1e-14
+    )
+
+    surface_markers = f.MeshFunction(
+        "size_t", mesh_fenics, mesh_fenics.topology().dim() - 1
+    )
+    surface_markers.set_all(0)
+    # Surface ids
+    left_id = 1
+    right_id = 2
+    left_surface.mark(surface_markers, left_id)
+    right_surface.mark(surface_markers, right_id)
+    ds = f.Measure("ds", domain=mesh_fenics, subdomain_data=surface_markers)
+
+    my_flux = SurfaceFluxSpherical("solute", right_id)
+    my_flux.D = f.Constant(2)
+    V = f.FunctionSpace(mesh_fenics, "P", 1)
+    expr = f.Expression(
+        f"{c_left} + {(c_right - c_left)/(r1-r0)} * x[0]",
+        degree=1,
+    )
+    my_flux.function = f.interpolate(expr, V)
+
+    my_flux.n = f.FacetNormal(mesh_fenics)
+    my_flux.ds = ds
+
+    computed_value = my_flux.compute()
+    # expected value is the integral of the flux over the surface
+    flux_value = float(my_flux.D) * (c_left - c_right) / (r1 - r0)
+    expected_value = -4 * math.pi * flux_value * r1**2
+    assert np.isclose(computed_value, expected_value)
+
+
+@pytest.mark.parametrize(
+    "azimuth_range", [(-1, np.pi), (0, 3 * np.pi), (-1, 3 * np.pi)]
+)
+def test_azimuthal_range_cylindrical(azimuth_range):
+    """
+    Tests that an error is raised when the azimuthal range is out of bounds
+    """
+    with pytest.raises(ValueError):
+        SurfaceFluxCylindrical("solute", 1, azimuth_range=azimuth_range)
+
+
+def test_soret_raises_error():
+    """
+    Tests that an error is raised when the Soret effect is used with SurfaceFluxCylindrical
+    """
+    my_flux = SurfaceFluxCylindrical("T", 1)
+    with pytest.raises(NotImplementedError):
+        my_flux.compute(soret=True)
+
+
+@pytest.mark.parametrize(
+    "azimuth_range", [(-1, np.pi), (0, 3 * np.pi), (-1, 3 * np.pi)]
+)
+def test_azimuthal_range_spherical(azimuth_range):
+    """
+    Tests that an error is raised when the azimuthal range is out of bounds
+    """
+    with pytest.raises(ValueError):
+        SurfaceFluxSpherical("solute", 1, azimuth_range=azimuth_range)
+
+
+@pytest.mark.parametrize(
+    "polar_range", [(0, 2 * np.pi), (-2 * np.pi, 0), (-2 * np.pi, 3 * np.pi)]
+)
+def test_polar_range_spherical(polar_range):
+    """
+    Tests that an error is raised when the polar range is out of bounds
+    """
+    with pytest.raises(ValueError):
+        SurfaceFluxSpherical("solute", 1, polar_range=polar_range)
+
+
+def test_soret_raises_error_spherical():
+    """
+    Tests that an error is raised when the Soret effect is used with SurfaceFluxSpherical
+    """
+    my_flux = SurfaceFluxSpherical("T", 1)
+    with pytest.raises(NotImplementedError):
+        my_flux.compute(soret=True)