diff --git a/stan/math/prim/fun/size.hpp b/stan/math/prim/fun/size.hpp
index a9005d97b39..34bb9f42358 100644
--- a/stan/math/prim/fun/size.hpp
+++ b/stan/math/prim/fun/size.hpp
@@ -34,6 +34,11 @@ inline int64_t size(const T& m) {
   return m.size();
 }
 
+template <typename T, require_tuple_t<T>* = nullptr>
+inline int64_t size(const T& /*t*/) {
+  return std::tuple_size<std::remove_reference_t<T>>{};
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/prim/fun/size_mvt.hpp b/stan/math/prim/fun/size_mvt.hpp
index 26b55e145d8..d766b500b9d 100644
--- a/stan/math/prim/fun/size_mvt.hpp
+++ b/stan/math/prim/fun/size_mvt.hpp
@@ -36,11 +36,21 @@ int64_t size_mvt(const MatrixT& /* unused */) {
   return 1;
 }
 
+template <typename TupleT, require_tuple_t<TupleT>* = nullptr>
+int64_t size_mvt(const TupleT& /* unused */) {
+  return 1;
+}
+
 template <typename MatrixT, require_matrix_t<MatrixT>* = nullptr>
 int64_t size_mvt(const std::vector<MatrixT>& x) {
   return x.size();
 }
 
+template <typename TupleT, require_tuple_t<TupleT>* = nullptr>
+int64_t size_mvt(const std::vector<TupleT>& x) {
+  return x.size();
+}
+
 template <typename StdVectorT, require_std_vector_t<StdVectorT>* = nullptr>
 int64_t size_mvt(const std::vector<StdVectorT>& x) {
   return x.size();
diff --git a/stan/math/prim/fun/vector_seq_view.hpp b/stan/math/prim/fun/vector_seq_view.hpp
index 8c76979927a..a28b916c7e3 100644
--- a/stan/math/prim/fun/vector_seq_view.hpp
+++ b/stan/math/prim/fun/vector_seq_view.hpp
@@ -8,12 +8,12 @@
 namespace stan {
 namespace internal {
 template <typename T>
-using is_matrix_or_std_vector
-    = math::disjunction<is_matrix<T>, is_std_vector<T>>;
+using is_matrix_or_std_vector_or_tuple
+    = math::disjunction<is_matrix<T>, is_std_vector<T>, math::is_tuple<T>>;
 
 template <typename T>
 using is_scalar_container = math::disjunction<
-    is_matrix<T>,
+    is_matrix<T>, math::is_tuple<T>,
     math::conjunction<is_std_vector<T>, is_stan_scalar<value_type_t<T>>>>;
 }  // namespace internal
 
@@ -76,7 +76,7 @@ class vector_seq_view<T, require_t<internal::is_scalar_container<T>>> {
  */
 template <typename T>
 class vector_seq_view<
-    T, require_std_vector_vt<internal::is_matrix_or_std_vector, T>> {
+    T, require_std_vector_vt<internal::is_matrix_or_std_vector_or_tuple, T>> {
  public:
   explicit vector_seq_view(const T& v) noexcept : v_(v) {}
   inline auto size() const noexcept { return v_.size(); }
diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp
index e04c7e87172..0bbc6b30391 100644
--- a/stan/math/prim/prob.hpp
+++ b/stan/math/prim/prob.hpp
@@ -111,6 +111,7 @@
 #include <stan/math/prim/prob/gamma_lcdf.hpp>
 #include <stan/math/prim/prob/gamma_lpdf.hpp>
 #include <stan/math/prim/prob/gamma_rng.hpp>
+#include <stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp>
 #include <stan/math/prim/prob/gaussian_dlm_obs_lpdf.hpp>
 #include <stan/math/prim/prob/gaussian_dlm_obs_rng.hpp>
 #include <stan/math/prim/prob/gumbel_ccdf_log.hpp>
diff --git a/stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp b/stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp
new file mode 100644
index 00000000000..2a068ced450
--- /dev/null
+++ b/stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp
@@ -0,0 +1,163 @@
+#ifndef STAN_MATH_PRIM_PROB_GAUSSIAN_COPULA_CHOLESKY_LPDF_HPP
+#define STAN_MATH_PRIM_PROB_GAUSSIAN_COPULA_CHOLESKY_LPDF_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/to_vector.hpp>
+#include <stan/math/prim/fun/rep_vector.hpp>
+#include <stan/math/prim/fun/vector_seq_view.hpp>
+#include <stan/math/prim/fun/size.hpp>
+#include <stan/math/prim/fun/size_mvt.hpp>
+#include <stan/math/prim/prob/multi_normal_cholesky_lpdf.hpp>
+#include <stan/math/prim/prob/std_normal_lpdf.hpp>
+#include <stan/math/prim/prob/std_normal_log_qf.hpp>
+
+namespace stan {
+namespace math {
+namespace internal {
+
+template <typename Ty_elem, typename Ttuple, std::size_t... I>
+auto apply_lcdfs_elem_impl(Ty_elem&& y_elem, Ttuple&& lcdf_tuple,
+                           std::index_sequence<I...>) {
+  auto&& lcdf_func = std::get<0>(lcdf_tuple);
+  return lcdf_func(std::forward<Ty_elem>(y_elem),
+                   std::get<I + 1>(std::forward<Ttuple>(lcdf_tuple))...);
+}
+
+template <typename Ty, typename Ttuple>
+auto apply_lcdfs(Ty&& y, Ttuple&& lcdf_tuple) {
+  return index_apply<std::tuple_size<std::remove_reference_t<Ttuple>>{}>(
+      [&y, &lcdf_tuple](auto... Is) {
+        return std::make_tuple(apply_lcdfs_elem_impl(
+            y[Is], std::get<Is>(lcdf_tuple),
+            std::make_index_sequence<std::tuple_size<std::remove_reference_t<
+                                         decltype(std::get<Is>(lcdf_tuple))>>{}
+                                     - 1>{})...);
+      });
+}
+}  // namespace internal
+
+/** \ingroup multivar_dists
+ * The log of the Gaussian copula density for the given y and
+ * a Cholesky factor L of the variance matrix.
+ *
+ * As the Gaussian copula requires inputs to be in the unit interval,
+ * users are required to provide a tuple (of the same size as y) where each
+ * element is a tuple containing a functor for calculating the LCDF and any
+ * additional parameters required by the LCDF functor.
+ *
+ * This version of the function is vectorized on y and the tuple of LCDF
+ * functor tuples.
+ *
+ * @param y A scalar vector
+ * @param lcdf_fun_tuple A tuple of tuples, where each inner tuple follows the
+ *   structure {functor, param1, param2, ...}, where the functor has signature
+ *   (y[i], param1, param2, ...) and returns the log CDF for the given y[i].
+ * @param chol The Cholesky decomposition of a variance matrix
+ * of the multivariate normal distribution
+ * @return The log of the gaussian copula density.
+ * @throw std::domain_error if LL' is not square, not symmetric,
+ * or not semi-positive definite.
+ * @tparam T_y Type of scalar.
+ * @tparam T_lcdf_fun_tuple Type of tuple of LCDF functor tuples.
+ * @tparam T_chol Type of cholesky factor.
+ */
+template <bool propto, typename T_y, typename T_lcdf_fun_tuple, typename T_chol>
+auto gaussian_copula_cholesky_lpdf(const T_y& y,
+                                   const T_lcdf_fun_tuple& lcdf_fun_tuple,
+                                   const T_chol chol) {
+  static constexpr const char* function = "gaussian_copula_cholesky_lpdf";
+
+  using T_y_ref = ref_type_t<T_y>;
+  using T_chol_ref = ref_type_t<T_chol>;
+  using T_lcdf_ref = ref_type_t<T_lcdf_fun_tuple>;
+
+  check_consistent_sizes_mvt(function, "y", y, "lcdf_fun_tuple",
+                             lcdf_fun_tuple);
+  const size_t size_mvt_y = size_mvt(y);
+  const size_t size_mvt_lcdf_tuple = size_mvt(lcdf_fun_tuple);
+  T_y_ref y_ref = y;
+  T_chol_ref chol_ref = chol;
+  T_lcdf_ref lcdf_tuple_ref = lcdf_fun_tuple;
+
+  vector_seq_view<T_y_ref> y_vec(y_ref);
+  vector_seq_view<T_lcdf_ref> lcdf_tuple_vec(lcdf_tuple_ref);
+  using T_return = return_type_t<
+      T_y, decltype(internal::apply_lcdfs(y_vec[0], lcdf_tuple_vec[0])),
+      T_chol>;
+
+  if (size_mvt_y == 0) {
+    return T_return(0);
+  }
+  const size_t size_vec = max_size_mvt(y, lcdf_fun_tuple);
+
+  const int size_y = math::size(y_vec[0]);
+  const int size_lcdf_tuple = math::size(lcdf_tuple_vec[0]);
+
+  // check size consistency of all random variables y
+  for (size_t i = 1; i < size_mvt_y; i++) {
+    check_size_match(function,
+                     "Size of one of the vectors of "
+                     "the random variable",
+                     math::size(y_vec[i]),
+                     "Size of the first vector of the "
+                     "random variable",
+                     size_y);
+  }
+  // check size consistency of all CDF tuples
+  for (size_t i = 1; i < size_mvt_lcdf_tuple; i++) {
+    check_size_match(function,
+                     "Size of one of the vectors of "
+                     "the LCDF functor tuple",
+                     math::size(lcdf_tuple_vec[i]),
+                     "Size of the first vector of the "
+                     "location variable",
+                     size_lcdf_tuple);
+  }
+
+  check_size_match(function, "Size of random variable", size_y,
+                   "size of LCDF functor tuple", size_lcdf_tuple);
+  check_size_match(function, "Size of random variable", size_y,
+                   "rows of covariance parameter", chol.rows());
+  check_size_match(function, "Size of random variable", size_y,
+                   "columns of covariance parameter", chol.cols());
+
+  for (size_t i = 0; i < size_vec; i++) {
+    check_not_nan(function, "Random variable", y_vec[i]);
+  }
+  check_cholesky_factor(function, "Cholesky decomposition of a variance matrix",
+                        chol_ref);
+
+  promote_scalar_t<T_return, std::vector<Eigen::VectorXd>> q(size_vec);
+  T_return lp(0);
+  for (size_t i = 0; i < size_vec; i++) {
+    const auto& y_i = y_vec[i];
+    const auto& func_i = lcdf_tuple_vec[i];
+
+    const auto& res = internal::apply_lcdfs(y_i, func_i);
+    const auto& u = index_apply<
+        std::tuple_size<std::remove_reference_t<decltype(res)>>{}>(
+        [&res, size_y](auto... Is) {
+          Eigen::Matrix<T_return, Eigen::Dynamic, 1> u_inner(size_y);
+          static_cast<void>(std::initializer_list<int>{
+              (u_inner[Is] = std::get<Is>(res), 0)...});
+          return u_inner;
+        });
+    check_bounded(function, "LCDF-transformed inputs", u, NEGATIVE_INFTY, 0);
+    q[i] = std_normal_log_qf(u);
+    lp -= std_normal_lpdf<propto>(q[i]);
+  }
+  const std::vector<Eigen::VectorXd> zero_vec(size_vec, rep_vector(0, size_y));
+  lp += multi_normal_cholesky_lpdf<propto>(q, zero_vec, chol_ref);
+  return lp;
+}
+
+template <typename T_y, typename T_lcdf_fun_tuple, typename T_chol>
+auto gaussian_copula_cholesky_lpdf(const T_y& y,
+                                   const T_lcdf_fun_tuple& lcdf_fun_tuple,
+                                   const T_chol chol) {
+  return gaussian_copula_cholesky_lpdf<false>(y, lcdf_fun_tuple, chol);
+}
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/test/unit/math/mix/prob/gaussian_copula_cholesky_test.cpp b/test/unit/math/mix/prob/gaussian_copula_cholesky_test.cpp
new file mode 100644
index 00000000000..60eb8f8bef8
--- /dev/null
+++ b/test/unit/math/mix/prob/gaussian_copula_cholesky_test.cpp
@@ -0,0 +1,31 @@
+#include <test/unit/math/test_ad.hpp>
+
+TEST_F(AgradRev, ProbDistributionsGaussCopulaCholesky) {
+  // Bind functors for compatibility with AD framework
+  auto f = [](const auto func1, const auto func2) {
+    return [=](const auto& y, const auto& args, const auto& sigma) {
+      auto lcdf_functors = std::make_tuple(
+          std::make_tuple(func1, args[0], args[1]), std::make_tuple(func2));
+      auto sigma_sym = stan::math::multiply(0.5, sigma + sigma.transpose());
+      auto L = stan::math::cholesky_decompose(sigma_sym);
+      return stan::math::gaussian_copula_cholesky_lpdf(y, lcdf_functors, L);
+    };
+  };
+
+  // y[0] ~ gamma(2, 1)
+  // y[1] ~ std_normal()
+  Eigen::VectorXd y1(2);
+  y1 << 1.0, 0.1;
+  Eigen::VectorXd args(2);
+  args << 2, 1;
+
+  auto func1 = [](const auto& y, const auto& shape, const auto& scale) {
+    return stan::math::gamma_lcdf(y, shape, scale);
+  };
+  auto func2 = [](const auto& y) { return stan::math::std_normal_lcdf(y); };
+
+  Eigen::MatrixXd Sigma22(2, 2);
+  Sigma22 << 2.0, 0.5, 0.5, 1.1;
+
+  stan::test::expect_ad(f(func1, func2), y1, args, Sigma22);
+}
diff --git a/test/unit/math/prim/prob/gaussian_copula_cholesky_lpdf_test.cpp b/test/unit/math/prim/prob/gaussian_copula_cholesky_lpdf_test.cpp
new file mode 100644
index 00000000000..b5a71abc6f2
--- /dev/null
+++ b/test/unit/math/prim/prob/gaussian_copula_cholesky_lpdf_test.cpp
@@ -0,0 +1,108 @@
+#include <stan/math/prim.hpp>
+#include <gtest/gtest.h>
+
+TEST(ProbDistributionsGaussianCopulaCholesky, EqualsMultiNormalCholesky) {
+  using stan::math::diag_pre_multiply;
+  using stan::math::gaussian_copula_cholesky_lpdf;
+  using stan::math::multi_normal_cholesky_lpdf;
+
+  Eigen::VectorXd y(2);
+  y << 1, 3;
+
+  Eigen::VectorXd mu(2);
+  mu << 0.1, 0;
+
+  Eigen::VectorXd sd(2);
+  sd << 2, 1;
+
+  Eigen::MatrixXd chol(2, 2);
+  chol << 1, 0, 0.5, 0.8660254037844386;
+
+  auto norm_lcdf = [](auto&& y, auto&& mu, auto&& sigma) {
+    return stan::math::normal_lcdf(y, mu, sigma);
+  };
+
+  auto std_norm_lcdf = [](auto&& y) { return stan::math::std_normal_lcdf(y); };
+
+  // y[0] ~ normal(0.1, 2)
+  // y[1] ~ normal(0, 1)
+  auto lcdf_functors = std::make_tuple(std::make_tuple(norm_lcdf, mu[0], sd[0]),
+                                       std::make_tuple(std_norm_lcdf));
+
+  double log_prob = stan::math::normal_lpdf(y[0], 0.1, 2)
+                    + stan::math::std_normal_lpdf(y[1])
+                    + gaussian_copula_cholesky_lpdf(y, lcdf_functors, chol);
+
+  double expected_log_prob
+      = multi_normal_cholesky_lpdf(y, mu, diag_pre_multiply(sd, chol));
+
+  EXPECT_FLOAT_EQ(log_prob, expected_log_prob);
+}
+
+TEST(ProbDistributionsGaussianCopulaCholesky, NonNormalMarginals) {
+  Eigen::VectorXd y(2);
+  y << 10, 6;
+
+  Eigen::MatrixXd chol(2, 2);
+  chol << 1, 0, 0.2, 0.9797958971;
+
+  auto gamma_lcdf = [](auto&& y, auto&& shape, auto&& scale) {
+    return stan::math::gamma_lcdf(y, shape, scale);
+  };
+
+  auto exp_lcdf = [](auto&& y, auto&& scale) {
+    return stan::math::exponential_lcdf(y, scale);
+  };
+
+  // y[0] ~ gamma(2, 1)
+  // y[1] ~ exponential(2)
+  auto lcdf_functors = std::make_tuple(std::make_tuple(gamma_lcdf, 2.0, 1.0),
+                                       std::make_tuple(exp_lcdf, 2.0));
+
+  double log_prob
+      = stan::math::gamma_lpdf(y[0], 2.0, 1.0)
+        + stan::math::exponential_lpdf(y[1], 2.0)
+        + stan::math::gaussian_copula_cholesky_lpdf(y, lcdf_functors, chol);
+
+  EXPECT_FLOAT_EQ(log_prob, -16.61005941);
+}
+
+TEST(ProbDistributionsGaussianCopulaCholesky, Errors) {
+  Eigen::VectorXd y(2);
+  y << 10, 6;
+
+  Eigen::VectorXd small_y(1);
+  small_y << 10;
+
+  Eigen::MatrixXd chol(2, 2);
+  chol << 1, 0, 0.2, 0.9797958971;
+
+  auto gamma_lcdf = [](auto&& y, auto&& shape, auto&& scale) {
+    return stan::math::gamma_lcdf(y, shape, scale);
+  };
+  auto exp_lcdf = [](auto&& y, auto&& scale) {
+    return stan::math::exponential_lcdf(y, scale);
+  };
+
+  auto lcdf_functors = std::make_tuple(std::make_tuple(gamma_lcdf, 2.0, 1.0),
+                                       std::make_tuple(exp_lcdf, 2.0));
+
+  auto small_lcdf_functors
+      = std::make_tuple(std::make_tuple(gamma_lcdf, 2.0, 1.0));
+
+  auto invalid_lcdf_functors
+      = std::make_tuple(std::make_tuple(gamma_lcdf, 2.0, 1.0),
+                        std::make_tuple([](auto&& y) { return y * 10; }));
+
+  EXPECT_THROW(
+      stan::math::gaussian_copula_cholesky_lpdf(y, small_lcdf_functors, chol),
+      std::invalid_argument);
+
+  EXPECT_THROW(
+      stan::math::gaussian_copula_cholesky_lpdf(small_y, lcdf_functors, chol),
+      std::invalid_argument);
+
+  EXPECT_THROW(
+      stan::math::gaussian_copula_cholesky_lpdf(y, invalid_lcdf_functors, chol),
+      std::domain_error);
+}