diff --git a/Migration.md b/Migration.md
index 305c099a1..3d73e92b6 100644
--- a/Migration.md
+++ b/Migration.md
@@ -5,6 +5,13 @@ Deprecated code produces compile-time warnings. These warning serve as
 notification to users that their code should be upgraded. The next major
 release will remove the deprecated code.
 
+## Gazebo Math 10.X
+
+### Breaking Changes
+
+1. **Color.hh**
+    + Fix: return type and behaviour change of member function `Color::operator[](const unsigned int _index)`. Now behaves like a mutator function, returning a mutable reference to the `Color` component (`float&`) instead of just the component value (`float`). In case of wrong index input, function returns a reference to `NAN_F`, and assigning any value to it has no effect. Refer to [#701](https://github.com/gazebosim/gz-math/pull/701) for further details.
+
 ## Gazebo Math 8.X to 9.X
 
 1. **SphericalCoordinates.hh**
diff --git a/include/gz/math/Color.hh b/include/gz/math/Color.hh
index b38beb8db..1bdd1b86b 100644
--- a/include/gz/math/Color.hh
+++ b/include/gz/math/Color.hh
@@ -168,7 +168,7 @@ namespace gz::math
     /// 3=alpha)
     /// \return r, g, b, or a when _index is 0, 1, 2 or 3. A NAN_F value is
     /// returned if the _index is invalid
-    public: float operator[](const unsigned int _index);
+    public: float &operator[](const unsigned int _index);
 
     /// \brief Array index operator, const version
     /// \param[in] _index Color component index(0=red, 1=green, 2=blue,
@@ -177,6 +177,14 @@ namespace gz::math
     /// returned if the _index is invalid
     public: float operator[](const unsigned int _index) const;
 
+    /// \brief Returns the Color component indexed by _index (accessible also
+    /// through const object)
+    /// \param[in] _index Color component index(0=red, 1=green, 2=blue,
+    /// 3=alpha)
+    /// \return r, g, b, or a when _index is 0, 1, 2 or 3. A NAN_F value is
+    /// returned if the _index is invalid
+    private: const float &GetRGBAfromIndex(const unsigned int _index) const;
+
     /// \brief Get as uint32 RGBA packed value
     /// \return the color
     public: RGBA AsRGBA() const;
diff --git a/src/Color.cc b/src/Color.cc
index 0759baa5c..8da83db2e 100644
--- a/src/Color.cc
+++ b/src/Color.cc
@@ -191,13 +191,21 @@ void Color::SetFromYUV(const float _y, const float _u, const float _v)
 }
 
 //////////////////////////////////////////////////
-float Color::operator[](const unsigned int _index)
+float &Color::operator[](const unsigned int _index)
 {
-  return (*static_cast<const Color *>(this))[_index];
+    return const_cast<float&>(
+            static_cast<const Color *>(this)->GetRGBAfromIndex(_index)
+            );
 }
 
 //////////////////////////////////////////////////
 float Color::operator[](const unsigned int _index) const
+{
+  return this->GetRGBAfromIndex(_index);
+}
+
+//////////////////////////////////////////////////
+const float &Color::GetRGBAfromIndex(const unsigned int _index) const
 {
   switch (_index)
   {
@@ -213,7 +221,13 @@ float Color::operator[](const unsigned int _index) const
       break;
   }
 
-  return NAN_F;
+  // This allows referencing a NAN_F value.
+  static float invalidParam;
+  // Prevents overwriting of invalidParam with non NAN_F
+  // value through wrong index.
+  invalidParam = NAN_F;
+
+  return invalidParam;
 }
 
 //////////////////////////////////////////////////
diff --git a/src/Color_TEST.cc b/src/Color_TEST.cc
index 348f8a60a..17c7e81c6 100644
--- a/src/Color_TEST.cc
+++ b/src/Color_TEST.cc
@@ -186,20 +186,36 @@ TEST(Color, Color)
   EXPECT_FLOAT_EQ(1.0f, clr.A());
 
   clr.SetFromHSV(300, 0.5, 1.0);
-  EXPECT_FLOAT_EQ(1.0f, clr[0]);
-  EXPECT_FLOAT_EQ(0.5f, clr[1]);
-  EXPECT_FLOAT_EQ(1.0f, clr[2]);
-  EXPECT_FLOAT_EQ(1.0f, clr[3]);
-  EXPECT_TRUE(std::isnan(clr[4]));
+  EXPECT_FLOAT_EQ(1.0f, clr.R());
+  EXPECT_FLOAT_EQ(0.5f, clr.G());
+  EXPECT_FLOAT_EQ(1.0f, clr.B());
+  EXPECT_FLOAT_EQ(1.0f, clr.A());
 
   clr.R() = 0.1f;
   clr.G() = 0.2f;
   clr.B() = 0.3f;
   clr.A() = 0.4f;
-  EXPECT_FLOAT_EQ(0.1f, clr[0]);
-  EXPECT_FLOAT_EQ(0.2f, clr[1]);
-  EXPECT_FLOAT_EQ(0.3f, clr[2]);
-  EXPECT_FLOAT_EQ(0.4f, clr[3]);
+  EXPECT_FLOAT_EQ(0.1f, clr.R());
+  EXPECT_FLOAT_EQ(0.2f, clr.G());
+  EXPECT_FLOAT_EQ(0.3f, clr.B());
+  EXPECT_FLOAT_EQ(0.4f, clr.A());
+
+  EXPECT_FLOAT_EQ(clr.R(), clr[0]);
+  EXPECT_FLOAT_EQ(clr.G(), clr[1]);
+  EXPECT_FLOAT_EQ(clr.B(), clr[2]);
+  EXPECT_FLOAT_EQ(clr.A(), clr[3]);
+  EXPECT_TRUE(std::isnan(clr[4]));
+
+  clr[0] = 0.5f;
+  clr[1] = 0.6f;
+  clr[2] = 0.7f;
+  clr[3] = 0.8f;
+  clr[4] = 0.9f;
+  EXPECT_FLOAT_EQ(0.5f, clr[0]);
+  EXPECT_FLOAT_EQ(0.6f, clr[1]);
+  EXPECT_FLOAT_EQ(0.7f, clr[2]);
+  EXPECT_FLOAT_EQ(0.8f, clr[3]);
+  EXPECT_TRUE(std::isnan(clr[4]));
 
   clr.Set(0.1f, 0.2f, 0.3f, 0.4f);
   clr = clr + 0.2f;
@@ -328,6 +344,12 @@ TEST(Color, ConstAndSet)
   EXPECT_FLOAT_EQ(clr.B(), 0.3f);
   EXPECT_FLOAT_EQ(clr.A(), 0.4f);
 
+  EXPECT_FLOAT_EQ(clr.R(), clr[0]);
+  EXPECT_FLOAT_EQ(clr.G(), clr[1]);
+  EXPECT_FLOAT_EQ(clr.B(), clr[2]);
+  EXPECT_FLOAT_EQ(clr.A(), clr[3]);
+  EXPECT_TRUE(std::isnan(clr[4]));
+
   math::Color clr2;
   clr2.R(0.4f);
   clr2.G(0.3f);
diff --git a/src/python_pybind11/src/Color.cc b/src/python_pybind11/src/Color.cc
index 9428fa26a..2dee0fd33 100644
--- a/src/python_pybind11/src/Color.cc
+++ b/src/python_pybind11/src/Color.cc
@@ -116,29 +116,33 @@ void defineMathColor(py::module &m, const std::string &typestr)
         &Class::SetFromABGR,
         "Set from uint32 ABGR packed value")
    .def("r",
-        py::overload_cast<>(&Class::R),
+        py::overload_cast<>(&Class::R, py::const_),
         "Get the red value")
    .def("g",
-        py::overload_cast<>(&Class::G),
+        py::overload_cast<>(&Class::G, py::const_),
         "Get the green value")
    .def("b",
-        py::overload_cast<>(&Class::B),
+        py::overload_cast<>(&Class::B, py::const_),
         "Get the blue value")
    .def("a",
-        py::overload_cast<>(&Class::A),
+        py::overload_cast<>(&Class::A, py::const_),
         "Get the alpha value")
    .def("r",
-        py::overload_cast<float>(&Class::R),
-        "Get the red value")
+        py::overload_cast<const float>(&Class::R),
+        "Set the red value",
+        "_r"_a)
    .def("g",
-        py::overload_cast<float>(&Class::G),
-        "Get the green value")
+        py::overload_cast<const float>(&Class::G),
+        "Set the green value",
+        "_g"_a)
    .def("b",
-        py::overload_cast<float>(&Class::B),
-        "Get the blue value")
+        py::overload_cast<const float>(&Class::B),
+        "Set the blue value",
+        "_b"_a)
    .def("a",
-        py::overload_cast<float>(&Class::A),
-        "Get the alpha value")
+        py::overload_cast<const float>(&Class::A),
+        "Set the alpha value",
+        "_a"_a)
    .def("__str__", toString)
    .def("__repr__", toString)
    .def("__copy__", [](const Class &self) {
@@ -148,7 +152,15 @@ void defineMathColor(py::module &m, const std::string &typestr)
      return Class(self);
    }, "memo"_a)
    .def("__getitem__",
-        py::overload_cast<const unsigned int>(&Class::operator[]));
+        py::overload_cast<const unsigned int>(&Class::operator[], py::const_),
+        "Array index operator, const version. Get the indexed component value",
+        "_index"_a)
+   .def("__setitem__",
+        [](Class &self, const unsigned int _index, const float value) {
+          self[_index] = value;
+        },
+        "Array index operator. Set the indexed component value",
+        "_index"_a, "value"_a);
 }
 }  // namespace python
 }  // namespace math
diff --git a/src/python_pybind11/test/Color_TEST.py b/src/python_pybind11/test/Color_TEST.py
index 0b2f09d0f..3afb3afd8 100644
--- a/src/python_pybind11/test/Color_TEST.py
+++ b/src/python_pybind11/test/Color_TEST.py
@@ -172,10 +172,26 @@ def test_color(self):
         self.assertAlmostEqual(1.0, clr.a())
 
         clr.set_from_hsv(300, 0.5, 1.0)
-        self.assertAlmostEqual(1.0, clr[0])
-        self.assertAlmostEqual(0.5, clr[1])
-        self.assertAlmostEqual(1.0, clr[2])
-        self.assertAlmostEqual(1.0, clr[3])
+        self.assertAlmostEqual(1.0, clr.r())
+        self.assertAlmostEqual(0.5, clr.g())
+        self.assertAlmostEqual(1.0, clr.b())
+        self.assertAlmostEqual(1.0, clr.a())
+
+        self.assertAlmostEqual(clr.r(), clr[0])
+        self.assertAlmostEqual(clr.g(), clr[1])
+        self.assertAlmostEqual(clr.b(), clr[2])
+        self.assertAlmostEqual(clr.a(), clr[3])
+        self.assertTrue(math.isnan(clr[4]))
+
+        clr[0] = 0.5
+        clr[1] = 0.6
+        clr[2] = 0.7
+        clr[3] = 0.8
+        clr[4] = 0.9
+        self.assertAlmostEqual(0.5, clr[0])
+        self.assertAlmostEqual(0.6, clr[1])
+        self.assertAlmostEqual(0.7, clr[2])
+        self.assertAlmostEqual(0.8, clr[3])
         self.assertTrue(math.isnan(clr[4]))
 
         clr.set(0.1, 0.2, 0.3, 0.4)