This paper introduces the asQ library for implementing ParaDiag parallel-in-time methods, focusing on applications related to weather and climate.  The paper is really well written, and I was impressed by both the asQ library design and numerical results.  The authors gave a nice intro on other relavent research with a comprehensive set of references.  They also did a good job of describing the asQ library and positioning it relative to existing methods and software.  The numerical examples are known to be difficult to solve parallel-in-time, yet excellent speedups were achieved in all cases.  The authors are aware of the current limitations of ParaDiag, but have several ideas for future research to expand its applicability.  I only have a few minor comments.

Minor comments:

- Page 2, paragraph before 1.1.  Change to "a survey of previous research".

- Page 6, after (13).  Change to "inverse x = P^{-1}b" or "inverse of Px = b"?  This is what the algorithms seems to be computing.

- Page 10, second sentence after (36).  Change to "backward Euler" (no "s").

- Page 10, "time_partition = [2, 2, 2, 2]".  This description grows with P_t.  This may not be a problem for a while (or ever, since roundoff grows with N_t), but I was curious.  Is there a P_t independent way to describe the partition in asQ?  Are there any other things about the interface that don't scale?

- Page 16, description of strong scaling.  I was really confused by this until I saw the subsequent comment about nonlinear problems and the "full timeseries of N_T timesteps" broken up into windows.  I'm not sure why you don't describe the linear problems in the same way, especially given that you use the term "windows" when discussing the numerical results.  Related to this, is the time given for the serial-in-time method in Figure 3 for 16,384 timesteps (which you could call N_T)?  Also in Figure 3, are you plotting T_p/N_t or (T_p/N_t)*N_w (where T_p is the time to solve N_t timesteps in a window)?  It looks like maybe you don't compute the full timeseries in the linear case, and instead compute only one window of size N_t and simply assume (reasonably) that the full N_T timesteps would be N_w times larger.  If this is all true, I think it would be less confusing to present it this way because it lines up exactly with the standard definition of strong scaling.  If I've misunderstood, then some additional detail is needed.

- Page 27, 5.2.1.  Delete one of "been" in the second sentence.

General comments.

The paper introduces asQ, a new software library built on top of the Firedrake software, for rapid testing and development of the ParaDiag parallel-in-time algorithm. While tailored towards earth system modeling, solved equations can be specified via the Uniform Form Language (UFL), making the framework quite general. The sections of the paper provide a detailed introduction into ParaDiag, the numerical algorithm, including a comprehensive review of the relevant literature. Then, the asQ library is described including the used space-time parallelization paradigm. The reader is walked through an example how to solve the heat equation with asQ to illustrate the steps that are needed. Wallclock scaling is shown for different numerical examples. While speedups for linear problems are extremely good, the fact that an averaged Jacobian needs to be used in ParaDiag results in much more modest speedups for nonlinear problems which are, however, very much in line with what other PinT methods deliver.

Overall, I think this is a strong, timely and well-written paper. Providing high-quality software libraries is exactly what is required to push parallel-in-time integration methods to broader use and there are still only few codes that do this. Therefore, asQ and the paper clearly fill an important need. Explanations in the paper are very clear and should be understandable for non-PinT-experts. The numerical results are convincingly linked to theory while the shown speedups are well explained using the provided performance model. I have only a few fairly minor issues that should be addressed before publication. The used code is provided open-source and numerical experiments can be reproduced via Python scripts and a complete installation is available as a Singularity container.

Specific comments.

1.  I am not entirely sure I understand why the transpose and thus all-to-all communication is needed. Is this a consequence of the way the method is implemented or is this intrinsic to ParaDiag?
   
   
   
   If I understand correctly, in Steps 1 and 3 on p. 6 only a (I)FFT in time is required for every spatial DoF, but the (I)FFT for each spatial DoF is independent from all others. Since the number of parallel copies of one spatial DoF (that is, N_t) is small(ish), is it not possible to store all the temporal copies of a given spatial DoF on the same node?
   
   
   
   Basically, let us say the spatial parallelization breaks down the spatial domain into four sub-domains Omega_1, Omega_2, Omega_3, Omega_4. Then say we use four parallel time steps. Cannot Node 1 hold the four temporal copies of Omega_1, Node 2 four copies of Omega_2 etc? This way, the (I)FFT in time can be computed without sending messages. This could even enable some "hybrid" space-time parallelization where the (I)FFT are parallelized with OpenMP or some other shared memory paradigm to avoid having to transfer full solutions in time. Is this possible in principle but not within the used framework, just very challenging to implement or is my understanding incorrect?
   
   
   
   I am not suggesting that the authors implement this but I would be interested in some more details on why the transpose is required and if, theoretically, there are ways around it.
2.  Eq. (24), I am not sure I understand what the Lip operator means. Is that the minimum L such that satisfies a Lipshitz condition for f(u) - Nabla f(u_tilde) * u ? But as what, as a function of u with fixed u_tilde? May be it is clearer to simply spell out the condition that the kappa needs to satisfy.
3.  Some of the figures are a bit hard to read in print. I would suggest slightly thicker lines and slightly larger markers. The yellow coloured lines are also sometimes difficult to see in print, may be some darker colour would help.
4.  There are a lot of parameters to keep track of (N_x, N_t, k_s, k_p, m_s, m_p, ...) and I found myself going back and forth a lot to find their definitions. Having them all summarized in a concise table for quick reference would help a lot.
5.  The description of the main steps of ParaDiag in the paragraph after Eq. (27) could probably be summarised in pseudo-code.

Technical corrections/minor comments.

1. The unit in which wallclock time is measured seems not to be stated and is missing from the axis labels in Figures.
2. There is a very recent preprint discussing integration of parallel-in-time functionality into the Nektar++ code. This effort probably deserves to be mentioned in the discussion of the (few) general purpose software packages supporting PinT.
   
   
   
   Xing, Jacques Y. and Moxey, David and Cantwell, Chris D., Enhancing the Nektar++ Spectral/Hp Element Framework for Parallel-in-Time Simulations. Available at SSRN: https://ssrn.com/abstract=5010907 or http://dx.doi.org/10.2139/ssrn.5010907